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INFORMAL FINAL REPORT ON 


BLADE FREQUENCY PROGRAM 
FOR NONUNIFORM HELICOPTER ROTORS, 

WITH AUTOMATED FREQUENCY SEARCH 
By S. Gene Sadler 

Rochester Applied Science Associates, Inc. 


SUMMARY 


A computer program has been implemented which determines the 
natural frequencies and normal modes of a lumped parameter model 
of a rotating, twisted beam, with nonuniforra mass and elastic 
properties. The end of the beam near the center of rotation may 
have one of four types of boundary conditions which are common to 
helicopter rotor systems; the outboard end has zero forces and 
moments, i.e., "free" boundary conditions. Six types of motion 
coupling may be modeled: fully coupled torsional-f latwise- 

edgewise motion, partially coupled torsional-f latwise motion or 
f latwise-edgewise motion, and uncoupled torsional motion, flatwise 
motion, or edgewise motion. Three frequency search methods have 
been implemented, including an automated search technique which 
allows the program to find up to the fifteen lowest natural fre- 
quencies without the necessity for input estimates of these fre- 
quencies by the user. 

This report contains computer program documentation informa- 
tion, and is intended to be sufficiently complete to be useful 
to both engineering users and programmers. 
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INTRODUCTION 


Normal modes and natural frequencies are used frequently in the 
analysis of elastomechanical systems for determination of instabili- 
ties and response to forcing functions. The use of these techniques 
for analysis of helicopter blades has been done for several years, 
with increasingly more complex models for both the helicopter rotor 
system and the aerodynamic load environment in which it operates. 

The increasing model complexities are tractible because of the im- 
plementation of simulation models on large-core high-speed digital 
computers. The need for more detailed analyses of air loads and 
blade response is particularly important for helicopters because of 
their normally complex aerodynamic environment, their elastomechani- 
cal characteristics, and because of increasingly more demanding per- 
formance requirements. Air loads and blade response are important 
from the viewpoint of acoustics, blade life, and vibration levels 
at the helicopter fuselage. Many helicopters have known (or anti- 
cipated) forcing function frequencies. In such cases, the blade 
frequencies and mode shapes themselves are of value in the blade 
design process. In other cases, the blade normal modes are of use 
as generalized coordinates in forced response calculations. The 
predicted response of normal modes is strongly dependent upon the 
predicted natural frequencies, and for systems such as helicopter 
blades, coupling between torsional and flatwise or flatwise and 
edgewise motions can cause significant shift in the predicted natu- 
ral frequencies, i.e., uncoupled or fully coupled frequencies. The 
most accurate (presumably the fully coupled) natural frequencies 
should be used in computing the forced response of such systems. 
However, there are cases where uncoupled or partially coupled modes 
and natural frequencies are also of interest. 

The theoretical methods and corresponding computer program 
which are discussed in this report provide an economical, easy-to- 
use, flexible tool for computing the normal modes and natural fre- 
quencies of lumped parameter models of a single helicopter blade (or 
other similar systems). It was expected that the user of this 
report a.nd program would have some familiarity with the methods of 
finding normal modes and natural frequencies of lumped parameter 
mass— elastic systems, and of their use in the design or analyses of 
such systems. However, the use of this report with its related 
references, together with elementary works on vibrating systems, 
should provide most potential users with sufficient background to 
make good use of this program. 
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THEORETICAL FORMULATION 


Matrix notation is common in many areas of application. The 
use of matrix methods in the field of linear elastomechanics has 
been highly developed. Transfer matrix methods provide a conven- 
ient formulation method and an efficient computational procedure for 
digital computer solution of physical problems which are described 
as a linear function of two variables, of which one is usually time 
and another is a spacial variable. Its application is especially 
well suited for determining the dynamic properties of linear elastic 
systems which can be considered to be a set of elements linked to- 
gether end to end, with little or no branching. Beams, shafts, and 
helicopter blades are examples of systems which are particularly 
well adapted to analysis by transfer matrix techniques. This tech- 
nique has been treated extensively in the literature, and the for- 
mulation, notation, and methods of solution discussed here are 
basically those of reference 1. The application of transfer matrix 
methods to determine the normal modes and natural frequencies of 
vibrating beams has been developed to the point that these charac- 
teristics may be rapidly and accurately determined for essentially 
any desired range of frequencies. Mode shape and frequency inaccur- 
acies, interdependence of one frequency on previously determined 
frequencies, and other difficulties which are often encountered in 
such calculations have been largely overcome by the methods dis- 
cussed herein. In particular, the basic method discussed herein is 
referred to as the modified transfer-matrix method, and its develop- 
ment and application are described in detail in Chapters 7 and 11 
of reference 1. 


Modified Transfer-Matrix Method 

Derivation of and examples of the use of general transfer 
matrices are discussed in detail in reference 1, particularly in 
Chapter 3. The method utilizes two relatively simple matrices, one 
which relates dynamical parameters at one point on a structure to 
those at another point on the structure, and one which is made up 
of those dynamical parameters of importance for the particular pro- 
blem which is under consideration; these are referred to as the 
transfer matrix, [T] , and the state vector, (z), respectively. The 
systems which are considered here have state vectors and transfer 
matrices as described in the following section. 

State vectors and transfer matrices .- The state vector at a 
point 5 is a vector, (z) ^ , whose elements are the deflections 

(linear and angular) and forces (shears and moments) which exist at 
that point. The state vector at another point, say j+1, is denoted 
by ( z )j +1 * The transfer matrix, [t] ^ , relates to state vector at 

the point j to that at the point j+1, and the relationship may be 
stated in the form 


3 



[T] j (z).. 


( 1 ) 


(z) 


j+1 


The state vector for fully coupled torsion-f latwise-edgewise motion 
of a helicopter blade has elements given by 


(z) - U,T,V,e,M z/ -Vy,-w,<J;,My,V z ) 


where v and w are deflections in the y and z directions, respective- 
ly/ <P , $ r and 8 are rotations about the x, y, and z axes, respec- 
tively, T is the torsional moment about the x axis, My and M z are 

bending moments about the y and z axes , respectively , and Vy and V z 

are shears in the y and z directions, respectively . The right- 
handed cartesian coordinate system and sign conventions for state 
vector elements are shown in figure 1 (taken from figure C-l of 
reference 1). The transfer matrices for many types of structural 
elements are given in reference 1 and other literature. Those which 
are discussed herein are repeated, for convenience, in APPENDIX B, 
along with a transfer matrix for a point (lumped) inertia and mass 
in a centrifugal force field. 
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Figure 1. 


Coordinate system and state vector sign 
convention 
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Blade elements and corresponding transfer matrices .- The trans- 
fer matrices, [T] , are normally a product of field or point matrices 
which correspond to idealizations of elements of the system which is 
being analyzed- The typical transfer matrix in this program is the 
product of the transfer matrices (which may be thought of as occurr- 
ing in order with increasing radius) for a rigid offset of the elas- 
tic axis in the chordwise direction, [C] , a point lumped mass and 
rotary inertia, [M] , a massless elastic length, [E] , and a rotation 
about the x-axis, [$] . (Detailed definitions of these matrices are 


given in Appendix B.) 


That is 



in the expression 
- IT] 1 (Z) D 


m. = in. [ej. [m] j [ c ]. 


( 2 ) 


Typical elements are sketched in figure 2 for the lumped para- 
meters model, and the signed elements are shown in element number 3 
with positive values. That is, e is positive for the mass center 
of gravity ahead of the elastic axis, h is positive for the elastic 
axis ahead of the radial coordinate line, and l is positive for an 

incremental change of elastic axis position which results in in- 
creased h with increased r. In the plan view of figure 2 blade 
twist is neglected, but when it exists the incremental change in 
twist, A<f>, is positive for increased nose-up twist with increased r. 



Figure 2. Typical blade elements with signed elements 
of lumped parameter ,model indicated 
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Increasingly more detailed models may be obtained by using more 
and shorter sections where properties change along the span. For 
lumped parameter models, it is suggested that at least two or three 
times as many degrees of freedom be allowed for a particular type 
of mode than the number of frequencies which are desired for that 
type of mode. For example, each element corresponds to one degree 
of freedom in torsion, so if four torsional normal modes and natural 
frequencies are desired, at least eight to twelve elements should be 
used in the model. Each element also corresponds to one flatwise 
degree of freedom and two edgewise degrees of freedom. The use of 
transfer matrices with distributed mass would reduce the number of 
elements needed for a given desired number of degrees of freedom. 
This would result in small savings in core storage and running time. 
The basis for and a discussion of some of the advantages of the use 
of such distributed property transfer matrices are contained in 
reference 2. 


Frequency determinant .- The application of transfer matrices in 
determining normal modes and natural frequencies of a rotating. blade 
is done through their use in defining the characteristic equation 
for the blade natural frequencies. Consider the blade shown in 
figure 3, which is made up of N elements, and for which the transfer 
matrices between adjacent state vectors are known. That is, 

(2) x = [T] x (Z) 0 , (Z) 2 = [T] 2 (Z) x , ..., (z) N = [T] n (z) n _ 1 , 

and the state vectors at the ends, ( z )g an( 3 (z) N are given for a 
cantilever - free beam, for example, as 

(z) Q = (0, T, 0, 0, M x , -V y , 0, 0, M , V z ) Q (3) 

(z) N = (<j> , 0, V, 0, 0, 0, -W, if,, 0, 0,) N (4) 



Figure 3. Blade model with N elements 
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In general, one-half of the elements of the state vectors at the 
ends of a structure are zero, and it is this fact that is used to 
obtain the determinental equation. Notice that the intermediate state 
vectors may be eliminated so that 

(z) N = IT] N ... IT] 2 IT] 1 (z) 0 (5) 

which may be written in the form 

(z) N = IP] Cz) q (6) 

where [P] is the product matrix, which results from the product of 
the transfer matrices. One of the advatanges of this method over 
some other methods is that this method can be implemented so that 
it uses less core storage than most other methods since only two 
matrices need to be in storage at one time, a matrix which is a 
product of two or more transfer matrices and the next transfer ma- 
trix to be used is defining the product matrix. 

Applying the boundary conditions, which are implicit in the 
state vectors at the ends of the beam, yields the determinant type 
of characteristic equation. That is, substituting equations (3) and 
and (4) into (6) results in 



and for a nontrivial solution 

D(«) = | Q | = 0 (7) 

where D(uj) is the determinant of the submatrix, [Q] , of the product 
matrix [P] , where [Q] is defined by 
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The rows of [Q] correspond to the zero elements in (Z)^ and the col- 
umns of [Q] correspond to the nonzero elements in (Z)q. This is the 


standard rule for defining IQ] from [P] . For a vibrating system, 
the transfer natrices, [T] , and therefore the product matrix, IP], 
the determinant matrix, IQ], and its determinant, D, are _ all func- 
tions of the vibration frequency, w. Values of to for which D(io) = 0 
are the natural frequencies of the blade. The definition of the 
determinant in the manner outlined above is straightforward but is 
often difficult to do with sufficient numerical accuracy that good 
mode shape quantities may be obtained. Important characteristics of 
the determinant as defined above are that D(co) is a polynomial in 
to 2 , has only positive roots, and is usually accurate numerically 
except very close to the natural frequencies of the system. 


Reduction of numerical accuracy problems .- A method which uti- 
lizes quantities that correspond to an approximate mode shape and 
mode shape correction terms in addition to an approximate frequency 
has been developed by Pestel and Leckie, and is explained in detail 
in reference 1, where it is referred to as the modified transfer- 
matrix method. In this method, an approximate mode shape, say U) Q / 

and correction terms, say [k] , are used in place of (Z)^; and the 

product matrix [P] , is replaced by an approximate mode shape and 
correction column matrix, [II], where the first column of In] Q is 

[ A ] Q , and the other columns of fn] Q are either zero or to, as dis- 
cussed in reference 1. That is, it is assumed (again using an N 
element cantilever - free beam as an example) that 
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and [H]q is made up of the column matrices which are used in the 
sum for (z)q. In general, (z ) is an (n, 1) matrix (or column) and 
tnj Q is an (n, n/2) matrix, where n is the number of physical quan- 
tities which are considered in the state vector. 

The transfer matrix process indicated in equation (5) is repeat- 
ed, but with ( z )q replaced by In] Q , and with the multiplications 

being done in the order indicated by the following steps 

[n], = IT]-, in]„ 

j- U 

[n] 2 = IT] 2 I nj 1 

• * 

• • 

In] 5 = IT] 5 [n] 4 

where it is understood that the columns in [n] i are to be added to 
obtain the column matrix (z) i « Since (z) N has zero elements, it is 
required that the [n]^ satisfy the appropriate set of homogeneous 
equations . 

For the specific example used in this discussion. 



where the columns which the k. multiply are the result of multiply- 
ing the corresponding columns in [n] Q by the [T] matrices. The 
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requirements that the [n ] N satisfy the boundary conditions as spe- 
cified by the zero elements in [ z ] . . is achieved if 
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It is shown in reference 1 that the A's in equation (10) are identi 
callv zero only if both the frequency, e , which was used in defining 
the It] and [ A] Q , which was part of [n ] Q , were a natural frequency, 

co n , and its corresponding state vector, I z ] ^ , respectively. (The 

A's in equation (10) are rarely exactly zero in practical cases.) 
Equation (10) may be thought of as an overdete mined set of equa- 
tions for the correction factors (k) . Quantities called remainders 
may be defined from these equations; they are more accurate than the 
determinant in the neighborhood of the natural frequencies, and 
their use is vital in overcoming most of the numerical accuracy pro- 
blems associated with use of the determinant type of characteristic 
equation which was discussed in the previous section. Their theo- 
retical derivation, examples of their use, and a procedure for their 
calculation are given in reference 1. The calculation procedure is 
outlined briefly here. (Chapter 7, section 3 of reference 1 also 
gives an example.) 

If the first row of equation (10) is eliminated, the (k) may 
be evaluated. That is 


Set b^ 
define 
from 


i — ! 


' K 51 
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(ID 


equal to this value of (b^ is called x in eq. 7-10 in ref. 
an alternate value for k^, say a^ (called y in reference 1) 


a l ^ K 22 k 2 + K 23 k 3 + K 24 k 4* ^ K 21 


( 12 ) 


1 ) 


Then define the remainder, , as 


R 1 ~ b l a l 


(13) 
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This remainder, and others defined by eliminating other rows of 
equation (10), etc., are the remainders which may be used in deter- 
mining the natural frequencies of the model. It is shown in refer- 
ence 1 that the remainders have zeros where the determinent has zeros 
and that the remainders, R^, are the sums of small numbers, a. and 

b^, i.e., that a^ and b^ have opposite signs. The determinant is 

usually the difference of large numbers. This fundamental differ- 
ence between the determinant and the remainders is the basis of the 
superior numerical accuracy of the modified transfer-matrix method. 
The remainders have apparent discontinuities, but this usually tends 
to be an inconvenience rather than a serious deficiency. Details of 
how the remainders and determinant may be used in various methods of 
determining the natural frequency of a system are discussed in a 
following section. Both the product matrix, [P] , and the state vec- 
tor and correction columns matrix, (n] , are computed and used by 
these methods. 


Mode shape definition .- As the determinant of the remainders 
approach zero, the approximate frequency and mode shape become more 
nearly correct, and at some point the frequency becomes sufficiently 
close to the numerical model's true natural frequency that it can be 
considered to be the natural frequency. The mode shape may be de- 
termined by using the corresponding values of the X's and k’s, and 
by adding the columns in the [n] , matrix to obtain (z)j. This 

method is the least subject to numerical inaccuracies. An alternate 
method uses the most accurate (z)q, and by successive multiplications 

by [T]^, [T] 2 / etc., defines (z)-^, (z) etc. The combination of 


all the (z) .'s contains the mode shape of each of the state vector 


quantities as a function of spanwise position. The zero mode shape 
quantities in (z)q are identically zero, but those quantities in the 

state vector at the other end of the structure, (z)^ which theoreti- 


cally should be zero are never identically zero in practice. Some 
indications of the accuracy of a calculated natural frequency and 
mode shape are the relative changes in magnitude in quantities which 
should be zero in the (z) N state vector to the same quantities at 


the next state vector position, (z) N _^. Normally, a (z) N "zero" type 

quantity should be at least three orders of magnitude smaller at the 
end than the same (z) N _^ quantity. Two types of mode shapes which 


are normalized to two different quantities are commonly used in the 
analysis of helicopter blades. One type of mode shape has a tip 
deflection quantity (subscript N) which has a value of unity. For 
shapes which include flatwise deflections, it is common to set v N 

equal to unity, otherwise ^ or (-w) N is usually set equal to unity. 


The second type of mode shape has a unity generalized mass, and is 
useful in programs which use normal modes and generalized coordi- 
nates to determine the response of the system to forcing functions. 
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Boundary Conditions 


As previously stated, the boundary conditions are inherently 
contained in the state vectors at the end of the structure. Excep- 
tions are structures which have restraints between the ends of the 
structure. These may be analyzed by this method also, and such 
systems are discussed in reference 1. There are boundary conditions 
which do not automatically result in one-half of the state vector 
elements being zero, but these may always be written in terms of 
another transfer matrix and a state vector which does have zero for 
one-half of its elements. 

Choice of various boundary conditions correspond to choice of 
various rows and columns of [P] to define [Q] . Of the many boundary 
condition combinations which are possible for the system of equation 
(5) , relatively few are of practical interest. The tip (subscript 
N) boundary conditions for a helicopter blade correspond to a force- 
free condition, or to a "free" end. The hub (subscript 0) boundary 
conditions usually involve hinged or clamped types of systems, and 
there are occasionally non-coincident hinges used on some helicopter 
rotor blades. Some rotors are continuous through the hub, e.g~. 
teetering rotors. These may be considered to be symmetrical beams , 
and one half of them analyzed, with both symmetrical and anti- 
symmetrical modes and frequencies required to describe the system. 
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Mode Coupling 


The previous discussion has dealt completely with fully coupled 
torsional-f latwise-edgewise normal modes and natural frequencies. 

It is sometimes helpful to use partially coupled or uncoupled modes 
in the analysis of helicopter blades and similar systems. Such data 
can be obtained by considering only a submatrix of each of the pro- 
duct matrices [P] and [n ] , the determinant matrix [Q] , and the state 
vector (z). The types of uncoupled modes which could be obtained 
from such manipulation of equations (5) through (10) , for example, 
would be uncoupled torsion, flatwise, and edgewise modes. Partially 
coupled modes would include torsional-f latwise, torsional-edgewise, 
and f latwise-edgewise. 

It should be noted that when uncoupled torsional modes are de- 
termined that no remainders are available, and that in general the 
number of remainders is given by (n/2)-l where n is the number of 
state vector elements. The state vector elements corresponding to 
the various types of mode coupling are given as follows 


Mode Types State Vector Quantities 

uncoupled torsional <f> , T 

v, 6, M , -V 

-w, ip, M y , V z 

<J>, T, v, 0, M z , -V 

v , 0 , M , -V , -w , ip , M , V 
z y r y z 

<j), T, V, 0, M , -V , -w, ip, M , V 

r, , , r z r y r r , y ' z 

It should be noted that partially coupled torsional-edgewise 
modes are not included as a mode coupling combination. Such modes 
cannot be determined by using a A<J> of 90° and switching values of 
Ely and EI z on input unless the mass eccentricity is nonzero in the 

y-direction. Such eccentricities are not included in the present 
model. Also, the neutral axis cannot be modeled to have a position 
other than in a plane normal to the rotor shaft which pases through 
the hub. That is, no precone is allowed in the blade model. 


uncoupled flatwise 
uncoupled edgewise 
coupled torsional-f latwise 
coupled f latwise-edgewise 
fully coupled 
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Frequency Search Methods 


Various types of frequency search methods have been used in 
connection with eigenvalue type problems in general. The three 
which are described here are some of several which have been used 
by the authors for searching for eigenvalues by hand or machine. 

They are not the best for all cases, but have proven to be useful 
in connection with the type of problem discussed here. One, which 
is applicable to elastomechanical systems with positive, real eigen- 
values is designed for use with a minimum of required user skill, 
computational time and cost, and is referred to as the automated 
frequency search method. The other methods are forerunners of this 
method, and are referred to as the frequency stepping and initial 
estimate search method. They have application to some systems for 
which the automated method should not be assumed to be reliable. 

Automated frequency search method .- The natural frequencies of 
the system are such that the numerically defined determinant, D(w) 
and the remainders, R^(to) , are zero when to is a natural frequency, 

say to . That is D(to )=0 and R. (to )=0. The numerical inaccuracy of 
n n in 

D(oj) may prevent it from being used to obtain a)^ with sufficient 

accuracy to allow definition of mode shapes; the k (« n ) are not so 

subject to numerical error, and have always been found to be satis- 
factory in determining to^ with sufficient accuracy to result in good 

mode shape determination. It can be shown that D(to) is a polynomial 

2 2 
in to , and has no negative (to ) roots. That is, 

2 22 22 2 2 2 
D(co)=K(w — w 2 ) ( w -w^)***^ "" ) • • • 


where K is some constant and the are natural frequencies, with 

all (oo^ 2 ) being positive. The R^(oj) are ratios of polynomials in 

u) 2 , and have discontinuities which correspond to roots of the deno- 
minator. These characteristics of D(w) and R^(u>) provide a basis 

for developing a completely automatic frequency search technique.^ 
That is, D(w) is a polynomial, and may be used in a Newton iteration 

2 

to obtain approximate values of w . No values of o> n are negative, 

so if the iteration is started near zero, the frequency which is 
found first should be the lowest, and should be approached from the 
low side, as is indicated in figure 4. Once the lowest frequency 

is determined, an auxiliary function, say F, (to) , may be defined 

i 2 2 

which has had that frequency removed. That is F^ ( w) =D ( w) / ( w -to^ ). 
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Figure 4. Iteration to lowest natural frequency 


The same process should yield the lowest frequency of F^(oj). The 

process may be continued for any desired number of frequencies. 
Numerical inaccuracy problems and dependence of one frequency (and 
its mode shape) upon previously determined frequencies may be elim- 
inated by using an appropriate remainder to do the final iterations 
for the . The remainders cannot be used during the early itera- 
tions because of their discontinuities. Care must be used in choos- 
ing one of the R^(w)'s which does not have a discontinuity near the 

desired w n . The choice of the best R. (to) is automatic, and depends 
upon the nearness of the co predicted by the R^(w)*s to that predict- 
ed by the D(to). Figure 5 shows the determinant, auxilliary function, 
and remainders for a typical case, for the lowest two natural 
frequencies. 
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(a) determinant, D, and auxiliary function, F 



(b) remainders 

Figure 5. Determinant, auxiliary function, and remainders 
for a typical automated frequency search. 


Frequency stepping .— Frequencies may be determined by this 

method by specifying a starting value, to^^, and a step size factor, 
f. The determinant and remainders are computed for each to, and to is 

incremented by defining to^^ = to^ + "^ (1 + f) until the determinant 

changes sign. The Newton iteration technique is then used with an 
appropriate remainder. After one natural frequency is determined, 
the frequency stepping continues, beginning with a frequency defined 
in terms of the frequency just determined, say to n , according to 


to 


(1) = 


to 


(1 + f) 



The method is illustrated in figure 6, and some specific short- 
comings are also indicated. While this option can give well-defined 
D(oo) and R^(w) curves, it may miss frequencies due to too large a 

step size (figure 6a) or too high an initial guess (figure 6b) , and 
may be inefficient if the step size is too small (figure 6c) . The 
step factor value, f, may be either positive or negative so the 
steps may yield increasing or decreasing w's. The danger of f being 
too large is more severe for large frequencies than for small, and 
the inefficiency of f being too small is worsened if the frequency 
spacing between natural frequencies is large. 


□ calculated determinant values 
O natural frequencies found jzr~X 


missed frequencies 


(a) step size too larae 


missed frequency 




(b) starting value too high 




(c) step size too small 


Figure 6. Frequency stepping option 






Frequency estimates method.— Frequencies may be determined by 
this method by defining estimated values of the natural frequencies 
and using a Newton iteration to search for the actual values, A 

second value of frequency to ^ , is defined (in order to use the 

Newton iteration) in terms of the estimated values, to ^ , as 


0 ) 


Cl) 



(1 + f) 


where ordinarily f is small. These values of w are used to start 
the Newton iteration. Where initial estimates are close to the 
actual model natural frequency it is possible to use the remainders 
for the iteration procedure. The primary shortcomings of this 
method are that the initial estimates may be incorrect, and that the 
remainder chosen to control the iteration may jump near one or more 
of the desired natural frequencies. Either of these shortcomings 
may result in duplication, or omission, of some frequencies. Typi- 
cal results of this method are illustrated in figure 7 . While the use 
of this method often requires exercise of skill and judgement by the 
user, it has consistantly been used to find natural frequencies of a 
number of helicopter rotor and propeller blades. It is less attrac- 
tive than the automated search option only because of the required 
user skills and the time spent in obtaining frequencies and mode 
shapes where the remainders all have discontinuities near a desired 
frequency . 



Figure 7. 


Initial estimates option 



COMPUTER PROGRAM IMPLEMENTATION 


A computer program has been developed which implements most of 
the features of the modified transfer-matrix method for determining 
frequencies and mode shapes of rotating, twisted, non-uniform beams. 
The program handles fully coupled, partially coupled, and uncoupled 
modes, four types of boundary conditions, and has three frequency 
search options. Details related to the formulation are contained in 
the section entitled THEORETICAL FORMULATION and in reference 1. 
Program input and output has been designed for a user who is some- 
what familiar with the problem of determining and utlizing normal 
modes and natural frequencies of lumped parameter models of physical 
systems. An effort has been made to simplify and minimize input and 
output for more efficient use in general, and for compatible opera- 
tion from computer terminals. The program consists of approximately 
1500 punched cards, executes (after compilation) in approximately 
33000 octal core storage locations, and requires approximately 0.25 
CPO seconds per iteration on a CDC 6600 computer when the program is 
compiled with an optimizing compiler. The number of iterations per 
frequency varies, but averages approximately 10. The program has 
accuracy limitations which render the presently implemented program 
useless for operation on computers that use less than 14 digit num- 
bers in floating point calculations. 

It has been anticipated that those portions of the program which 
a user may be apt to modify, or may want to understand to a greater 
degree than others include those portions of the program related to 
(1) the definition of boundary conditions, (2) the definition of the 
type of mode coupling, and (3) the implementation of the frequency 
search methods. These areas of probable interest are discussed 
briefly. It should be expected that a large portion of currently 
existing and future helicopter rotor systems will be able to be 
analyzed in terms of single blade natural frequency and mode shape 
characteristics without modifying the existing program. These brief 
discussions are intended to help the user better understand the cap- 
abilities of the program, and should be useful in the event that 
program modifications are desired. 

The functions of the various program routines are stated briefly, 
with references to detailed discussions given where they exist. A 
computer program symbol dictionary and a table of dimensioned vari- 
ables which defines how dimensions may be changed for different 
computer simulation models are presented. Finally, a description of 
computer program input and output are given. Details related to the 
punched mode shape output are given since this is expected to be one 
area where users may desire to make changes. 
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Boundary Condition Options 


The program, as presently implemented, has four types of hub 
boundary condition options available. These are defined in sub- 
routines BN DRY 1 or BNDRYX, and are used as specified bv the input 
parameter, IBC. 

For the value of IBC indicated they are: 

(1) edgewise hinge outboard of flatwise hinge (fully articu- 
lated rotor) , 

(2) clamped flatwise and edgewise (cantilever hub) 

(3) pinned flatwise and clamped edgewise, and 

(4) clamped flatwise and pinned edgewise. 

Of these, (3) and (4) are used to calculate the antisymmetrical and 
the symmetrical modes, respectively, of a teetering rotor. The 
cantilever boundary conditions may be used with inboard element 
stiffness adjusted so as to represent flexures, etc., if desired. 

The input number, IBC, or the mode shape quantities which are print- 
ed at the inboard end of the model may be used to determine the type 
of boundary conditions which exist. 

Torsional deflections have been assumed to be clamped at the 
hub, with any control system flexibility lumped with the appropriate 
section's torsional stiffness. The effective section stiffness, 

EIX , is defined in terms of the section length, £ , the control 

spring stiffness k, and the actual section stiffness, EIX^, as 

EIX = EIX / (1 + EIX / ilk), 
e a a 

The inboard blade element has its inboard end located a distance 

X , from the center of rotation, and X , = 0 is allowed. The 
root root 

edgewise hinge of the fully articulated rotor is located at the out- 
board end of section 1, and that section may have zero length. The 
outboard state vector has zero force type quantites, i.e., the tip 
is "free" and these are set automatically by the program. 

The boundary conditions subroutines, BNDRY1 and BNDRYX de- 
fine the hub state and the vectors ( z ) 0 an< 3 (n) Q with one quantity, 

the hub torque, set equal to unity during the iterations. When a 
natural frequency has been obtained, the hub state vector is multi- 
plied by a factor to yield a tip state vector with one quantity set 
equal to unity, or with a mode shape generalized mass of unity, as 
was discussed previously. The boundary conditions at the tip are 
defined in agreement with the input values of NTYPE , and always 
correspond to a force-free condition at the tip. 
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The boundary conditions at the hub as defined m the hub state 
vector, ( z ) c , for the different boundary conditions which are pre- 
sently implemented are given by the following table, (Note that the 
hub has been assumed to be clamped in torsion in all cases! . The 
pinned flatwise-pinned edgewise boundary conditions are for the flat 
wise hinge at the inboard end of the first section and the edgewise 
hinge at the outboard end of the first section. All other boundary 
conditions are applied at the inboard end of the first section. 


TABLE 1. 

STATE VECTOR ELEMENTS AT HUB FOR VARIOUS 
BOUNDARY CONDITIONS 


IBC Values 

1 

flatwise restraint 
edgewise restraint 

pinned 

pinned 



Hub State Vector Quantities 


algebraic names 


elements used by computer program 


4 

0 

0 

0 

0 

T 

1.0 

1.0 

1.0 

1.0 

V 

0 

0 

0 

0 

0 

X 1 

0 

X 1 

0 

M z 

0 

X 1 

0 

X 1 

-V 

y 

X 2 

X 2 

X 2 

X 2 

-w 

0 

0 

0 

0 

Tp 

X 3 

0 

0 

X 3 

M 

y 

0 

X 3 

X 3 

0 

V 

z 

X 4 

X 4 

X 4 

X 4 














Mode Coupling Options 


The computer program will treat several models, including sub- 
sets of the fully coupled torsional-f latwise-edgewise modes. The 
partially coupled and uncoupled modes are obtained by specifying the 
appropriate input, NTYPE; the modes and frequencies are computed by 
using matrix partitioning to obtain the desired data. Uncoupled 
torsional, flatwise, or edgewise modes and frequencies may be cal- 
culated, as may partially coupled torsional-f latwise or f latwise- 
edgewise modes and frequencies. 

The hub and tip boundary conditions for the type of coupling 
defined by NTYPE are automatically chosen correctly for the boundary 
conditions of the type which are specified by IBC. The mode shape 
quantities which are not part of the partially coupled or uncoupled 
modes are printed or punched as zeros. The values of NTYPE, and 
resulting type of mode are 


NTYPE 


type of mode 


1 

2 

3 

4 

5 

6 


uncoupled torsional 

uncoupled flatwise 

uncoupled edgewise 

coupled torsional-f latwise 

coupled f latwise-edgewise 

coupled torsional-f latwise-edgewise 


Program control variables for definition of the proper matrix parti- 
tioning and tip boundary conditions, are done in the main line of 
each value of NTYPE. A default for NTYPE > 6 or NTYPE < 1 is imple- 
mented, and gives NTYPE = 6 results as the default value. No error 
message is printed in this case. 

It should be noted that partially coupled torsional-f latwise 
modes are not included as a possible type of mode coupling. The 
user should be warned against attempting to obtain such modes by 
making a 90° rotation about the x-axis by inputting a A$ of 90 with 
appropriate changes in EIY and EIZ, etc., since the mass eccentri- 
city in the flatwise direction and the rotary mass moment of inter- 
tia about a chord line are assumed to be zero. As a result, the 
normal type of mass coupling can not be modeled for partially 
coupled torsional-edgewise modes. 
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Frequency Search Methods 


Three frequency search options are available with the present 
program. The input variable NFIND determines the option used. In- 
put of NFIND = 1, 0, or -1 results in use of the automated search 
method, the estimates method or the stepping method, respectively. 

Of these methods, the automated search technique is expected to be 
the most useful. All input variables must have some numerical val- 
ues. However, all those variables defined by NAMELIST are not used 
for each of the frequency search options. Those not used for a 
particular option are given below, and required or suggested values 
for those which are used are given in the section on Description of 
Input. 

Automated frequency search .- The automated frequency search 
begins with two internally defined initial guesses near u> =1.0, and 
uses a Newton's iteration to determine the lowest natural frequency, 
An associated function, F^( u ), is then defined in the form 

F ^(u)) = D(oj) / (to 2 - to x 2 ) . 

Thus, the zero of D(to) at to^ has effectively been removed from F^(w). 
Two values of D(to) which were determined in the search for io^ were 

2 2 

saved or re-computed and are divided by the appropriate (to - co^ ) 

factor. These provide, by extrapolation, an initial estimate for 
the next higher natural frequency, to^, which is again found by a 

Newton's iteration procedure. 

In the search for the third natural frequency, F.(io) is defined 
as 


F 2 (to) = D (to ) / [(to 2 - to 1 2 ) (to 2 - to 2 2 )] 


and extensions to higher frequencies are obvious. The process is 
continued until the desired number of frequencies are determined or 
until all frequencies below an input cutoff frequency are determined, 
along with their corresponding mode shapes. This describes the over- 
all approach. In the iteration to find a zero of D(uj) , both the 
determinant and the associated function, F(to), may become inaccurate, 
but the remainders, R^(w), will not. Also, use of F^(oj) to evaluate 

frequencies above the first would mean that such frequencies would 
depend on all previous frequencies and numerical inaccuracies. There- 
for, once the iteration using F(w) closely approximates the pro- 
gram chooses an appropriate R^(u), and concludes the iteration using 
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that R. ( to ) . The nearness of to to a natural frecruencv is measureid bv 
1 . th 

the relative change in frequency from one iteration, say the 1 , to 

the next. If | ( to (i+1 > - co (i) ) 2 |/ uo (l+1) - _ . 316228 NEXP ° , a remain- 

der is chosen and used for the subsequent iterations. The R^(io) s, 

of course, are idependent of previous w^'s, and are numerically 

accurate, but they cannot be used during the entire process since 
they are not polynomials in to, and display apparent discontinuities. 
The remainder which is chosen for the subsequent iterations is that 
which predicts the next frequency value closest to that predicted bv 
the determinant, using remainder and determinant values computed 
during the previous two iterations. 


The iteration sequence for a particular natural frequency 
will be stopped if the number of iterations, NTIME, exceeds NMAX. 

In this event, the natural frequency counter NFREQ will be set 
equal to NUMF, and the program will go to the next model. (This 
is necessary to avoid incorrect definition of the auxiliary function 
and subsequent erratic behavior of the automated frequency search.) 
Therefore, it is suggested that NMAX be set to at least 30 for use 
of the automated frequency search option. 


Input variables not used when NFIND = +1 are NOOR, PRCNT, and, 
if NPDF-0, OMINT. 

Frequency estimates .- Starting estimates of natural frequencies 
are input to the program in the OMINT array. A second frequency 
value factor, PRCNT, for use in starting the Newton iteration, is 
also input, and should be small, about 1.0 or less. This option is 
implemented so that one remainder, as specified by the input vari- 
able NOOR, is used for the entire frequency range. The values of 
OMINT need not be in any order, and may have any numerical value. 
Input variables not used when NFIND = 0 are NPDF and OMLIMT. 

Frequency stepping search.- A step factor, PRCNT, and frequency 

2 

starting value, OMINT (1), are input and the values of to are stepped 
according to 

(co (l+1) ) 2 = (<jo (l) ) 2 (1 + PRCNT / 100) 

until the determinant changes sign. The program then begins a Newton 
iteration, using the remainder specified by the input index, NOOR. 

The initial value for the next frequency is defined from the natural 
frequency according to 


(1 + PRCNT / 100) 



The process is repeated until the specified number of frequencies 
are found. The value of PRCNT to be used for a particular case 
depends on the intended use of the option. If a scan over a large 
frequency range is desired, PRCNT is large. If a detailed scan over 
a limited frequency range is desired, then PRCNT is small. Input 
variables not used when NFIND = -1 are NPDF, OMINT(I) for I>1, and 
OMLIMT . 
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Program Routine Functions 


The overall flow of the program is indicated by figure 8. The 
program BLADE reads and writes input; sets program control con- 
stants; calls subroutine BNDRY1 or BNDRYX to establish hub boundary 
condition; performs transfer matrix multiplications; defines deter- 
minants and remainders (using subroutine MATINV to evaluate determi- 
nants and solve simultaneous sets of equations for the definition 
of remainders) , does the frequency search options of frequency esti- 
mates and frequency steppina (using subroutine INTERP) , and calls 
subroutine FINDM for the automatic search; calls subroutine SHAPE to 
define mode shapes when a natural frequency is found or when the 
limit on the number of iterations is reached; and controls the num- 
ber of frequencies per model and the number of models per run. 

Subroutine SNTFR sets program control equivalent to that given 
by NTYPE = 6 for input values of NTYPE which are not in the accept- 
ed range of 1 thru 6. 

Subroutine BNDRY1 defines the hub boundary conditions for IBC=1. 
This option gives modes for a fully articulated rotor with the edge- 
wise hinge at the outboard end of section 1 and with the flatwise 
hinge at the inboard end of section 1. Subroutine BNDRYX defines 
the hub boundary conditions for IBC = 2, 3, and 4 respectively. 

IBC = (2) gives clamped-clamped modes, IBC = (3) gives pinned- 
flatwise clamped-edgewise modes; and IBC = (4) gives clamped- 
flatwise pinned-edgewise modes. 

Subroutine TEMAT defines the transfer matrix for a typical 
blade element; see equation (2) and Appendix B for the definition 
of this matrix. Note that for various values of NTYPE , the blade 
element is assumed to be rigid (have infinite stiffness) in torsion, 
edgewise bending, or flatwise bending (as given in the following 
table) if NX, NY or NZ , respectively, are unity; otherwise, the 
blade has flexibility. These parameters make the calculation some- 
what more efficient for partially coupled or uncoupled modes. 


NTYPE 

NX 

NY 

NZ 

type of coupling 

1 

0 

1 

1 

uncoupled torsion 

2 

1 

0 

1 

uncoupled flatwise 

3 

1 

1 

0 

uncoupled edgewise 

4 

0 

0 

1 

coupled torsion-f latwise 

5 

1 

0 

0 

coupled f latwise-edgewise 

6 

0 

0 

0 

coupled torsion— f latwise-edgewise 

Values of 

NX, NY, 

or 

NZ of unity 

cause special printed messages 

after the 

run title 

is printed. 



BLADE 


BLADE 


BNDRY1 

or 

BNDRYX 

TEMAT 


BLADE 


MATINV 


FINDM 

and/or 

INTERP 

FINDM 

or 

BLADE 


BLADE 



Figure 8. Program Flow 
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Subroutine MATINV solves sets of simultaneous equations and 
estimates the determinant of the coefficient matrix. See reference 
3 for more details. 

Subroutine FINDM controls the automatic frequency search and 
printout of data during the search. The determinant value D(1.0) is 
stored during the first calculation for a particular model and is 
not recomputed. The next frequency value is always 1.001, and this 
value and the corresponding values of the determinant, auxiliary 
function, and remainders are printed as iteration number 1. Subse- 
quent frequency estimates and corresponding data are printed until 
the frequencies of successive iterations satisfy the relationship. 


(w 


(i) 2 


- 0 ) 


(i-1) 2 


) / « 


(i-1) 2 < 


- .316228 


NEXPO 


The remainder of the iterations are done by using the "best" avail- 
able remainder value. (Except in the case of uncoupled torsional 
modes, for which there are no remainders and the determinant retains 
good accuracy. ) The "best" remainder is that which perdicts the 
next frequency estimate closest to that predicted by the determinant, 
using the two previous frequency estimates. The remainder chosen is 
identified on the output, and subsequent iterations use only that 
remainder. Values of the frequency determinant and remainders but 
not the auxiliary function are printed during the subsequent itera- 
tions. Usually only one or two iterations are necessary to satisfy 
the frequency convergence criterion as given by 

(J 1 * 2 - a* 1 - 1 * 2 ) / oV' 1 ’ 2 i 0.1 NEXPO 


Upon satisfaction of this criterion, subroutine FINDM resets program 
control variables to start searching for the next highest frequency 

and stores co 1 in the OMINT array. 

When one or more frequencies are stored in the OMINT array, 
they are used to define the auxiliary function as described in the 
THEORETICAL FORMULATION section. (The numerical values stored in 
OMINT are actually the squares of the frequencies.) 

Subroutine INTERP does a linear extrapolation (or interpolation) 
for the Newton iteration technique. 

Subroutine SHAPE controls the calculation, printing, and punch- 
ing of the mode shapes. (Subroutine SHAPE calls subroutine BNDRY1 or 
BN DRY X and TEMAT. ) 

Standard library functions or subroutines which are used include 
DAYTIM, EXP, SINCOS, SQRT , and TANH. 
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Computer Program Symbol Dictionary 

Computer program symbols are given here in alphabetical order 
with brief definitions or descriptions. Program routine and func- 
tion names are described in the previous part of this section. 
Symbols described as temporary storage quantities may be used for 
such purposes in more than one location within a program routine and 
within more than one routine. Details of their use may be obtained 
from the program listing if necessary. When temporary storage 
quantities are used in only one routine, that routine is indicated. 
Input quantities are preceeded by an asterisk. The user may refer 
to the Description of Input portion of this section for details con- 
cerning input quantity physical dimensions and allowable range of 
values, where appropriate. 


Computer 

Algebraic 

Definition or Description 

name 

symbol 


A 

[tt] 

mode shape + correction matrix and 
temporary storage array in SHAPE 

AMAX 


temporary storage quantity in MATINV 

ASAV 


temporary storage for [tt] 

AVMAS 

m 

lumped mass 

A107 

M 10 , 7 

element of [M] , see APPENDIX B 

Aij 

M. . 

elements of [M] , see APPENDIX B 

B 


temporary storage vector 

BAND 


temporary storage array in BNDRY1 

BK 


temporary storage quantity in BNDRYX 

BND 


arqument equivalent of PRD, in BNDRY1 
and BNDRYX 

BPL 


arqument equivalent of A, in BNDRY1 
and BNDRYX 

B2 

a l 

remainder quantity, see equation (12) 

C 


temporary storage quantity 

CAP SI 

Y 

angle between blade element's chord 
line and plane perpendicular to rotor 
shaft 

CCPS 

COS Y 


CL 


temporary storage quantity in TEMAT 

CL2 and CL3 


temporary storage quantities in TEMAT 

CNVRT 


tt/18 0 
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Computer 

Algebraic 

Definition or Description 

name 

symbol 


*CPOMG 

ft 

rotor speed 

CPSQ 

ft 2 


CSSI 

temporary storage for CCPS in TEMAT 

CSSQ 

(ft cos T ) 2 


Cl thru C8 


temporary storage quantities in TEMAT 

C101 


temporary storage quantities for C in 
BLADE 

DATE 


run date, defined by calling routine 
DAYTIM 

DET 

D(to) 

determinant 

DETERM 


argument in MATINV, matrix determinant 

DET1 

D(1.0) 

determinant 

DT 

|K| 

determinant of [K] of equation (10) 

DTOLD 


storage for determinant from previous 
frequency iteration 

DTSAV 


temporary storage for DET in FINDM 

El 

El 

temporary storage for EIX, EIY , or EIZ 
in TEMAT 

*EIX 

El 

X 

torsional rigidity 

*EIY 

El 

v 

chordwise bending stiffness 

*EIZ 

El 

z 

flatwise bending stiffness 

*EL 

l 

length 

*ELC 

l 

c 

control link length 

*ELZ 

A£ 

z 

elastic axis chordwise offset 

*EMAS 

m 

mass 

*EPS 

£ 

mass eccentricity 


or e 

£ i = (£ i + £ i-l )/2 

EPSEM 

em 


EPSN 

e N 

tip section eccentricity 

EPST 

temporary storage for EPS in TEMAT 

EPSX 


temporary storage for EPS in BLADE 

ERR 

[k] 

kappa vector, see equation (10) 

ERR101 


temporary storage for ERR in BLADE 
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Computer 

name 

FACTR 

F10 

F2 

F6 

GAMMA 

GNMS 


Algebraic 

symbol 


(1+f) 0 ' 5 


Y 

2 

Y 


H 


h 


HI 

I 

*IBC 

ICOLUMN 

ID 

ID1 

II 
IKJ 
IM1 
INDEX 
*INSEC 
I ROW 
*ITAB 
IX 


Definition or Description 

(1 + PRCNT / 1000) 
mass load factor 
inertia load factor 
mass load factor 
Y = * 0? T / El) 0 * 5 

mode shape generalized mass 


offset of elastic axis from radial 
coordinate 

temporary storage for H in TEMAT 


matrix index 

boundary condition control 
matrix index 

argument in MATINV, not used 

argument in MATINV, but used 

matrix index 

matrix index 

matrix index 

matrix index 

number of blade sections 

matrix index 

blade property table input control 
matrix index 


J 

JK 

JKII 

JKL 

JX 


matrix index 
matrix index 
matrix index 
matrix index 
matrix index 


K 

KDl 

KD2 


matrix index 

argument in BNDRYX, boundary condition 
control 

boundary condition control 
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Computer 

name 

KD3 

KD4 

L 

LL 

LI 

M 

Ml 

N 

NBND 

NC 

NCI 

NC2 

NC3 

NC4 

NC5 

*NEXPO 

NF 

*NFIND 

NFLAG 

NFLG 

NFREQ 

NGYRO 

NIBC 

NINT 

*NMAX 

NN 

*NOOR 

NOUT 

*NPCH 

*NPDF 

NPIB 

NR 


Algebraic Definition or Description 

symbol 

boundary condition control 
boundary condition control 

matrix index 
matrix index 
matrix index 

matrix index 

argument in MATINV, not used 

argument in MATINV, not used 
tip boundary condition indices 
number of columns in [T] , NC2-NC1-1 
(NR1 + l)/2 
NR2/2 
NC2 - NCI 
NC- 1 
NC- 2 

convergence criterion limit 
matrix index 

frequency search option selector 
frequency search control 
frequency search control 
q frequency counter 

= 0 

boundary condition control 
matrix index 

maximum number of frequency iterations 
internal program control 
remainder number used for iterations 
output tape number, = 6 
punched output control 
number of known frequencies 
argument in BNDRY1, program control 
NR2 - NR1 + 1 
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Computer 

name 

Algebraic 

symbol 

Definition or Description 

NR1 


index of first row in [T] 

NR2 


index of last row in [T] 

NR3 


NR-1 

NR4 


NR- 2 

NSEC 


sINSEC 

NSM1 


NSEC-1 

NSP1 


NSEC+1 

NSTEP 


frequency search control 

NTAPE 


frequency search control 

NTIME 


iteration counter 

NTRL 


frequency search control 

*NTYPE 


type of mode coupling 

*NUMF 

Snax 

number of frequencies to be found 

NX 


torsional rigidity control 

NXl 


argument in MATINV, not used 

NY 


chordwise bending rigidity 

NZ 


flatwise bending rigidity 

N1 


argument in MATINV, matrix size 

OM 


2 . 

next estimate for to in INTERP 

OMD 


auxiliary function factor 

OMEG 


temporary storage for OMEGA 

OMEGA 

(0 

natural frequency 

*OMINT 


frequency array 

*OMLIMT 

OMLSS 

OMLST 

2 (l) 
V i+1) 

03 

2 • 

limit on to magnitude 

2 

current iteration value for w 

2 

previous iteration value for to 

OMNXT 


temporary storage for OMOR in FINDM 

OMOLD 


auxiliary function factor 

OMOR 


temporary storage for OMSQ in FINDM 

OMP 


OMD* OMOLD 

OMQ 


temporary storage for OMRD 

OMRD 


magnitude of difference between values 
2 

of to predicted by D and various R's 
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Computer 

Algebraic 

Definition or Description 

name 

symbol 


OMSAV 


temporary storage for OMSQ in FINDM 

OMSQ 

2 

to 


*PHI 

A $ 

blade element twist 

PHIX 


storage for PHI 

PI 

IT 

3.1415926 

PIVOT 


diagonal matrix element in MATINV 

*PRCNT 


frequency search step size factor 

PRD 

[P] 

product matrix 

PTX 


axial centrifugal force in element 

PTZ 

P X 

T z 

chordwise centrifugal force in element 

RNEW 

R (i) 

remainders for current iteration 

ROLD 

rCx-D 

remainder for previous iteration 

S 

sinh y/y anc ^3 
(sinh y~y)/y 


SAV 


temporary storage for [P] 

SAVXX 


temporary storage for SAV 

SCSI 

sin'Fcos'y 


*SIZER 

0 

o 

collective pitch at blade root 

SIZERX 


storage for SIZER 

SL 

i /El 


SL2 

a 2 /ei 


SL3 

£ 3 /EI 


SNSI 

sin ¥ 


SNSQ 

( Q sin ¥) 2 


SOC 


variable used in determination of 
sinh y 

SUMAS 

m 

in . = (m . + m. . ) /2 
i l i-l 

SV1010 


temporary storage for SAV in BLADE 

SWAP 


temporary matrix element storage 

T 

[T] 

transfer matrix 

*THETC 

6 

c 

pitch link control angle 

THETCX 


storage for THETC 
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Computer 

name 

Algebraic 

symbol 

Definition or Description 

*TITLE 


run title 

*TMSL 


torsional mode punched output control 

TMR 

1 Z (1)/Z (3) 1 N 

torsional mode limit ratio 

VNORM 


mode normalization factor 

X 

r 

radial coordinate 

VT 

aj. 


temporary storage for X 

*XINR 

I 

x 

torsional inertia 

*XROOT 

x root 

length from shaft center to first 
flexible blade element 

YA 


storage for D or R 

YB 


storage for D or R 

*YINR 

I 

y 

chordwise inertia 

Z 

(z) 

state vector 

ZERMO 


temporary storage for Z(9)^ in BNDRY1 

ZSV 


temporary storage for Z in SHAPE 
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Dimensioned Variables 


Program variable dimensions may be grouped into three groups: 
those which depend on the maximum number of elements in the model, 
MS , those which depend on the maximum number of natural frequencies 
MF, and those which have a constant dimension size. Program vari- 
ables and the appropriate dimensions are given for each routine, in 
alphabetical order, with dimensions indicated as MS, MF , or the 
fixed number necessary. 


Routine Variable name (dimension) 

BLADE A (10,5) 

ASAV (10,5) 

AVMAS (MS ) 

B (10) 

C (10) 

CAPSI (MS) 

CCPS (MS) 

C101 (10,1) 

EIX (MS) 

EIY (MS) 

EIZ (MS) 

EL (MS ) 

ELZ (MS) 

EMAS (MS) 

EPS (MS) 

ERR (10) 

ERR101 (10,1) 

NBND (5) 

OMINT (MF) 

PHI (MS) 

PHIX (MS) 

PRD (10,5) 

RNEW ( 4 ) 

ROLD (4) 

SAV (10,5) 

SV1010 (10,10) 

T (10,10) 

TITLE (19) 



Routine 


Variable 


XINR (MS 
YINR (MS) 


same as BLADE except 

Delete: A 

ASAV 

C 

NBND 

OMINT 

PHIX 

PRD 

RNEW 

ROLD 

SAV 

T 

Add: BAND (10,5) 

BND (10,5) 
BPL (10,5) 

SAV (10,7) 


FINDM 

OM (4) 
OMINT (MF) 
OMRD (4) 


RNEW (4) 
ROLD (4) 

INTERP 

none 

MATINV 

A (10,10) 


B (10,1) 


INDEX (10,3) 


BNDRY1 

and 

BNDRYXj 



Routine 


Variable name (dimension) 
SHAPE same as BLADE except 

Delete: ASAV 


ERR 


NBND 


OMINT 


PHIX 

KNEW 

ROLD 

SAV 


TITLE 


Add: Z (10) 

ZSV (10) 


NBND (5) 


same as SHAPE except 
Delete: A 


C 


Description of Input 


Input parameters are provided to the computer on punched cards. 

A description of the input is given in this section in tabular form, 
and includes the card number, input format (and FORTRAN statement 
number) , and the computer program variable name, type, definition, 
and permissible range of values and/or physical dimensions where 
applicable. Physical dimensions are given and output labeled m 
English units. Any other system of units is acceptable as long as 
blade physical properties are consistent. However, output will 
remain labeled in English units unless changed. 

The FORTRAN Statement numbers and formats used to read input 
parameters are as follows: 


SNTFR 

TEMAT 
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401 FORMAT (211, 12, 19A4) 

402 FORMAT (I2,F6.5,9E8.0) 

and other parameters are read using NAMELIST, where variables given 
in the NAMELIST "IN" are punched in any order or format on cards 
which begin with a "$IN", which end with a "$" punched after all 
data, which have no punches in the first column of any card, and 
with variables separated by commas. 

Permissible values or ranges of values for input parameters are 
given where essential for program operation; other remarks concern- 
ing program and model limitations as discussed in the section on 
PROGRAM USE should also be taken into consideration. All input 
variables which begin with letters I through N are integer numbers , 
and are read without decimals and must be right justified in the 
input field. All other input variables are floating point variables, 
and are read in F or E format as described, except TITLE, which may 
contain alphanumeric data, and uses A format. The arrangements of 
input parameter description is as follows. 


(Card No.); Format No., Parameters defined 

Column Computer Algebraic Permissible values, physical 

name symbol dimensions, etc. 

(1) ; 401, Program control and model title 

1 IBC Boundary conditions control 

=1, pinned-pinned , with lead-lag 
hinge at outboard end of 
first element and flap hinge 
at inboard end of first 
element; 

=2, clamped flatwise and edge- 
wise , 

=3, pinned flatwise and clamped 
edgewise, or 

=4, clamped flatwise and pinned 
edgewise. 
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2 


I TAB 


Blade property table input con- 
trol 

/0, reads blade properties, or 
=0, uses previous model's blade 
properties. Mote that ITAB=0 
results in printed output values 
of PHI and EPS which correspond 
to the PHI in radians and the 
averaged values of EPS, e, rather 
than the "input" numbers as are 
printed for ITAB=1 


3,4 

INSEC 

N 

Number of elements in the blade 
model (INSEC 30 stops execution) , 
and this may be used as a "stop 
card" in the data deck. 

5-80 

TITLE 


Title of run (alphanumeric data) 

( 2 thru 

INSEC + 1) 

? 402, Blade 

element lumped parameter values 

1-2 

J 

j 

Blade element lumped number, in- 
creases with increasing radial 
coordinate 

3-8 

EL 

a 

Blade element length, feet 

9-16 

EIX 

EI 

X 

Element torsional stiffness, 
lb-ft 2 

17-24 

eiy 

EI 

y 

Element edgewise stiffness, 
lb-ft 2 

Column 

Computer 

Algebraic 

Permissible values, physical 


name 

symbol 

dimensions, etc. 

25-32 

EIZ 

EI 

z 

Element flatwise stiffness, 
lb-ft 2 

33-40 

XINR 

I 

X 

Element torsional inertia about 
2 

c. g. , lb-sec -ft 

41-48 

YINR 

I 

y 

Element chordwise inertia about 
c.g. (=XINR for models where 
numerical values are unavailable) 

49-56 

EMAS 

m 

2 

Element mass, lb-sec /ft 
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57-64 

PHI 

A <j> 

65-72 

EPS 

e 

73-80 

ELZ' 

l 

z 


(INSEC + 2) ; NAMELIST/IN/, Program 

(Order and columnar spacing is not 
numerical value must be included i: 


Element twist, degrees 

Mass eccentricity, positive for 
mass center of gravity forward 
of elastic axis, ft. 

Change in offset of elastic axis 
from a radial line drawn through 
the center of rotation, positive 
for average elastic axis location 
of present element ahead of that 
for previous element, feet 

control and operating conditions 

defined; the name and at least one 
l the NAMELIST/IN/ data) . 


NTYPE Type of mode coupling; 

=1, uncoupled torsional; 

=2, uncoupled flatwise; 

=3, uncoupled edgewise; 

=4, partially coupled torsional- 
f latwise; 

=5, partially coupled flatwise- 
edgewise; 

= 6, fully coupled torsional- 
f latwise-edgewise 

NMAX i m Maximum number of iterations 

allowed during the search for 
one natural frequency value, 
suggested values: 14 for NFIND^l, 
30 for NFIND=1. 


Computer Algebraic 
name symbol 

NEXPO 


NUMF 


q. 


max 


Permissible values, physical 
dimensions, etc. 

Convergence criteria limit para- 
meter, suggested value: 10 to 12 

Limit on the number of frequen- 
cies found for the current blade 
model, NUMF<15 due to dimension 
sizes (MF parameter in Dimension- 
ed Variables portion of the 
COMPUTER PROGRAM IMPLEMENTATION 
section) 
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NOOR 

NPCH 

CPOMG 

SIZER 

THETC 

ELC 

XROOT 

Computer 

name 

PRCNT 


NFIND 


\ Remainder number to be used for 

NFIND= - l,or 0; Range: 0<NOOR£ 

C (number of elements in state 
vector) /2-1) . Values of NOOR 
out of this range are superceded 
by the maximum value allowed for 
the current mode types. 

Punched output control: 

=0, no punched output 
>0, punched output: unity general- 
ized mass type of mode shape 
and frequency data, using 
formats 916 and 930 of sub- 
routine SHAPE. 

fj Rotor speed, rad/sec. Values of 

CPOMG £ 0 are allowed. 

¥ Blade root collective pitch 

° angle, degrees, any value 

allowed 

0 pitch link control rod angle 

c with plane perpendicular to 

rotor shaft, degrees. If no 
value is known, a suggested 
value is 0.0. 

1 Control link length, ft 
c 

x , Length from center of rotation 

to inboard end of first blade 
element, XROOT > 0 may be used. 


Algebraic 
symbol 


Permissible values, physical 
dimensions, etc. 


Step factor in frequency search 
methods, suggested values: 
not used for NFIND=1; 

£0.01 for NFIND=0 ; 

1 £ PRCNT £ 30 for NFIND=1 , 
depending on known or expected 
frequency spacing. See figure 
6. If PRCNT=0 is input, it is 
revised to PRCNT-10. 

Frequency search control 
=1, automated search option 
=0, initial estimates option 
=~i, frequency stcppino option 



NPDF Number of previously determined 

frequencies to be eliminated 
during automated frequency search 
(NFIND=1) . These values of co 

n 

are stored in the OMINT array. 


OMINT 


TMSL (V Z 3 > n 

max 


OMLIMT a) 2 

max 


Natural frequency estimates array 


NFIND 

No. of estimates 

-1 

1 

0 

greater of NUMF and 

+1 

greater of NPDF and 


Torsional mode size limit. If 
TMSL- | (Z(l)/Z(3j | the punched 
output is controlled by NPCH, 
otherwise no punched output will 
be obtained. Suggested values: 
TMSL > 1000 for no limit on 
punched modes 

TMSL = 1.0 for elimination of 
punched torsional 
modes 

Limit on the size of to , beyond 

which no natural frequencies are 
desired. Program stops after 

the first value of OMEGA 2 which 
exceeds OMLIMT. 


Alphabetic list of input parameters .- An alphabetic list of the 
computer program input variable names together with their respective 
data card group and definition are given here for convenience. The 
card group is one of the following: (1) title, (2) blade element 

properties, or (3) NAMELIST /IN/. 


Name 


CPOMG 

EIX 

EIY 

EIZ 

EL 

ELC 

ELZ 

EMAS 

EPS 


Card Definition 

group 


3 rotor speed 

2 torsional stiffness 

2 edgewise stiffness 

2 flatwise stiffness 

2 length 

3 control link length 

2 offset of elastic axis from radial 

coordinate 
2 mass 

2 mass eccentricity 
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IBC 
INSEC 
I TAB 


1 boundary condition control 

1 number of elements in blade model 

1 blade property input control 


Name 

J 


Card Definition 

group 

2 blade element number 


NEXPO 

NFIND 

NMAX 

NOOR 

NPCH 

NPDF 

NTYPE 

NUMF 


3 

3 

3 

3 

3 

3 

3 

3 


convergence criterion limit parameter 
frequency search control 
maximum number of iterations per 
frequency 

remainder number control 
punched output control 
number of previously determined 
frequencies 

type of mode coupling control 
limit on the number of frequencies 


OMINT 

OMLIMT 


3 

3 


natural frequency estimates array 

2 

limit on size of u 

n 


PHI 

PRCNT 


2 twist increment 

3 frequency search step factor 


SIZER 


3 


blade root collective pitch angle 


THETC 

TITLE 

TMSL 


3 pitch link angle 

1 run title (alphanumerical) 

3 torsional mode limit switch 


XINR 

XROOT 


2 torsional inertia 

3 length from center of rotation to 
first blade element 


YINR 


2 


chordwise inertia 
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Description of Output 


Output is printed in blocks for each model. Punched output of 
unity generalized mass type of mode shapes results for NPCH 0, and 
is not obtained for NPCH = 0. (No control over printed output is 
implemented.) The output for each model is arranged in three groups: 
(1) program and blade identification, including run title, run date, 
blade operating conditions and blade properties; (2) program control 
parameters; and (3) frequency search data and resulting mode shapes. 
All input parameters are written out, and labels and arrangement are 
intended to maximize clarity and convenience to the user. 

The program and blade model group of input parameters has the 
arrangement shown in APPENDIX D, page 107. The only variations from 
that format are those due to special messages concerning rigid models, 
as controlled by NTYPE, and those noted on page 40 when ITAB=0. 

The program control group of input parameters has values of 
OMINT printed only when they are used, otherwise that output is 
invarient. Typical output is shown in APPENDIX D on page 111. 

The frequency search data and mode shape parameters will be 
repeated until either NUMF or OMLIMT controls cause the program to 
stop. This output group normally requires, for each frequency, 
approximately one page for frequency search data and one or two pages 
for mode shape data. The data for both the frequency searches and 
mode shapes is given in columns, with the columns labeled, except 
that the' natural frequency in radians per second and as a multiple of 
rotor speed and the generalized mass of the mode are printed across 
the page immediately below each mode shape data set. 

The punched output results only if NPCH 0. The punched out- 
put corresponds to the printed mode shape data for the mode with 
unity generalized mass except that: (1) the data at station number 

"0" is not punched, and (2) a frequency in radians per second, the 
first 64 alphanumerical characters of 'the run title, and the number 
of the frequency, as given by the frequency counting index, NFREQ, is 
punched on a card preceding the mode shape data. Mode shape data is 
punched with two cards for each station. The first card contains 
the state vector elements z(3), z(7), z(l), z(8), and z(4) in that 
order in format 5(E14.7, IX) and the second card contains z(2), z(5), 
z(6), z(9), and z (10) in that order in format 5E15.7. Page 120 of 
APPENDIX D is a listing of cards corresponding to the printed input 
given on page 103 of APPENDIX D. The punched data may be used with 
no changes in the program discussed in reference 5. It is expected 
that other programs would require different formats, but such 
changes should be straightforward. All such changes would involve 
only two existing output statements in subroutine SHAPE. One is 
located three cards after FORTRAN statement number 17 and punches 
the natural frequency, title, and frequency number, and the other is 
given by FORTRAN statement number 51 and the following continuation 
card, and punches the state vector quantities. 


45 



PROGRAM USE 


This program was designed to compute the fully coupled normal 
modes and natural frequencies of helicopter blades, but can be used 
to analyze other similar physical systems or subsystems. The rota- 
tional speed may be positive or zero. The program is applicable to 
beams with an elastic axis which runs approximately parallel to a 
radial coordinate. The beams may have principal elastic axes which 
are twisted, may have nonuniform mass, inertia, and stiffness prop- 
erties and noncoincident mass center and elastic axes. The neutral 
axes for torsion and flexure should be approximately coincident. 
Steady pitch and twist angles are modeled, but steady coning and 
lagging angles are not modeled. A variety of hub end boundary condi- 
tions are available, but the tip end is assumed to be "free", i.e. 
force type quantities are zero and displacement type quantities are 
nonzero at the tip. The types of motion coupling which may be 
modeled are: coupled torsional-f latwise-edgewise motion, partially 
coupled torsional-f latwise motions and f latwise-edgewise motions, 
and uncoupled torsional, flatwise, and edgewise motions. 

The model used for the beam is a lumped parameter model which 
has lumped point mass and inertia and uniform massless elastic sec- 
tions. Twist of the beam is modeled by finite incremental twists at 
the mass location. The offset of the elastic axis from a radial 
line is also modeled by finite incremental offsets at the mass loca- 
tion. Centrifugal forces in the offsets are included (i.e. large 
steady forces in the chordwise direction) , but the offset lengths 
are assumed to be rigid. The elastic axes of the flexible sections 
are modeled to be parallel to the radial coordinate. Blade precone 
and prelag are not included in the model. 

This type of model does not provide a consistent lower or upper 
bound to the theoretical continuous model natural frequency, but 
should provide acceptable values for use in blade design and forced 
response analysis. A discussion of the errors in frequency estimates 
which are introduced by lumped parameter models is presented in 
reference 2. 

Rotating or non-rotating beams with the outboard end free, heli- 
copter rotor blades, propellers, fan and turbine blades are examples 
of the types of systems which can be modeled by this program. 
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Development of Blade Property Input Data 


Blade properties are usually available in graphical or tabular 
form, and lumped parameter models may be developed from such data in 
the following manner. The section lengths are chosen so as to per- 
mit uniform element subsections to model nonuniform blade properties. 
Usually the largest spanwise variations occur in mass or stiffness 
properties. It is suggested that a reasonable minimum number of ele- 
ments for a uniform blade is ten, but more elements may be needed if 
more than three of four flatwise modes and frequencies are desired. 
The blade mass eccentricity, if not available, may be assumed to be 
zero. The blade element mass, m, inertia, I and I , and incremental 

twist, A$; are the average values per unit length times the section 
length. The torsional and bending stiffness El , El , and El , 

« ^ y ^ 

the elastic offset from the previous section l , and the eccentri- 

z 

city, e, are the average values for the element, subject to the 
following general rules, which have been found to be useful. 


(1) Localized regions of relatively low torsional or 
bending stiffness, such as that indicated in 
figure 9 (a) , should be modeled as accurately as 
possible, especially if they occur near the center 
of rotation. 

(2) Localized regions of relatively large stiffness 
• may be neglected. 

(3) Localized variations in mass and inertia properties, 
eccentricity, and elastic axis location may be 
averaged, but should not be ignored. Large, highly 
localized mass distribution, such as indicated in 
figure 9(b) should be accurately modeled. 

(4) If a choice must be made, it is usually better to 
have inboard elements shorter than outboard elements . 

(5) Chordwise inertia, I , are often not given, and in 

such cases they should be made equal to the torsional 

inertia, I . 

x 
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Boundary Condition Options 


The hub boundary conditions should be chosen to most closely 
approximate the actual system. It should be noted that the fully 
articulated rotor must have its edgewise hinge "outboard" of the 
flatwise hinge, even though the length, £., between them may be zero. 
That is, the flatwise hinge is fixed to the hub and the edgewise 
hinge "flaps". The location of the torsional pitch bearings is 
arbitrary, and is determined by the section which has it EI^. value 

adjusted to account for a point torsional spring. If spanwise loca- 
tion is important, this section length may be made very small. 


Pinned boundary conditions result in exactly zero moments at 
the pin. Bending flexures may be modeled by clamped boundary condi- 
tions with appropriate section bending stiffnesses. Other types of 
boundary conditions may be modeled in a similar manner or new bound- 
ary condition options added in subroutine form. 


Most of the complications of the BNDRYl and BNDRYX subroutines 
arise because of the necessity to be able to model six combinations 
of motion coupling. It is suggested that a user may implement 
different boundary conditions, other than those presently imple- 
mented, most easily in the following manner. 

(1) Define the boundary conditions in terms of zero and nonzero 
state vector elements. 

(2) Define an intermediate transfer matrix composed of zeros 
and ones which when multiplied by one of the present bound- 
ary conditions gives the desired zero and nonzero state 
vector elements, for fully coupled motion. 

(3) Make program changes which perform the necessary multipli- 
cation of the existing PRD and A matrices referred to in 
(2) (which are provided by the BNDRYl or BNDRYX subrou- 
tines) prior to their multiplication by the standard transfer 
matrix, T, (which is provided by subroutine TEMAT ) . 


Mode Coupling Options 

Fully coupled modes and natural frequencies are preferred over 
partially coupled or uncoupled modes and natural frequencies for the 
analysis of blade response. This is primarily because the blade 
response is strongly dependent upon natural frequency values, and 
natural frequencies of fully coupled modes are generally more accurate 
than natural frequencies of partially coupled or uncoupled modes. 
However, even the fully coupled modes and frequencies do not account 
for mechanical coupling between blades except for teetering rotors, 
and no torsional coupling between blades is modeled for any boundary 
conditions. Mechanical coupling between blades and nonuniform 
swashplate stiffness can be important factors, and can be analyzed 
by the method discussed in reference 4. 
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Frequency Search Options 


Suggested uses of the frequency search options are discussed 
here, with a brief summary of their (intended) operating character- 
istics. Specific details of program input for these options are 
given in the section on Input. 

Automated frequency search.- Ordinarily the lower natural fre- 
quencies of a helicopter rotor blade (and many other systems) are of 
more interest than the higher natural frequencies. When all (or a 
majority) of the lower natural frequencies (up to a maximum of 15) 
are of interest, the use of the automated frequency search option 
provides a method for their determination which required a minimum 
of user skill in developing program control input. Because of its 
features, it would usually be the most economical option with regard 
to both labor and computer costs. 

In the event that some natural frequencies of a system are known 
and it is desired that these not be calculated, they are input in 
the OMINT array, with the value of NPDF set equal to the number of 
these previously determined frequencies. The program treats these 
as roots of D(u), and neither these frequencies nor their correspond- 
ing mode shapes are determined. (This feature permits the user to 
obtain, in a "second" run, some of the lowest fifteen natural fre- 
quencies which may have been missed in a "first" run due to time 
limits or which may have had poor mode shapes due to improper con- 
vergence limits, etc.) 

One of the limitations of this option is that the method as 
implemented does not allow the user to obtain frequencies above the 
lowest 15. However, this could be resolved by increasing the dimen- 
sions of those variables which depend on the maximum number of fre- 
quencies allowed. 

Initial estimate option .- When frequencies are not too closely 
spaced, or where approximate values are known, initial estimates 
provide a convenient means of obtaining their values. While this 
option has been somewhat obsoleted by the automated search option, 
it may be utilized in connection with that option in determining 
frequencies which the automated option missed due to the limit on 
the number of frequencies, or which did not yield good mode shapes 
due to convergence limits. Also, frequencies higher than the lowest 
15 may be obtained by using this option. 

Frequency stepping option.- This option may be used to compute 

determinants and remainders as a function of to , or to obtain de- 
tailed information over a limited frequency range. It is largely, 
obsolete due to the automated frequency option, but has been retain- 
ed for possible user needs. 
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Use of Output 


The output from the program consists of printed and punched 
data. The printed data includes all input quantities. These quan- 
tities are useful in checking input data during calculations, and in 
properly identifying and using previously generated output. The 
above quantities are labeled and should require no user modifications 
unless program input is changed. During the iteration process, 

* values of the iteration counter, current frequency estimate, deter- 

minant, auxiliary functions, and remainders are printed. Once a 
natural frequency is obtained or the limit on the number of itera- 
tions is reached, a mode shape is printed. The mode shape accuracy 
can be evaluated from the relative size of the force-type quantities 
at the tip as compared to those quantities at the next inboard state 
vector position. 

As was previously discussed, two types of mode shapes are print- 
ed out. The first has a normalized tip deflection quantity of unity 
and the second has a unity generalized mass. Where normal modes and 
natural frequencies are to be used to compute the forced response of 
blades in programs such as that discussed in reference 5, it is con- 
venient to have this second type of mode. The program has the option 
to punch the frequency and mode shape for NPCH = 1 , for use in such 
blade response calculations. The mode shape quantities at the in- 
board end of the first section are not included in the punched out- 
out, and are not normally used in calculations of blade response. 

In many design processes, natural frequencies are adjusted by 
modification of blade properties so as to avoid known forcing func- 
% tion frequencies. Mode shapes and natural frequencies in their 

printed output form would be useful for such design processes. 

The variables read from the title card and NAMELIST /IN/ include 
most of those which are expected to vary during blade analysis and 
design processes and for use in forced response calculations. If a 
user desires to redefine some blade properties, such a property 
could be changed readily by reading it in the NAMELIST /IN/ group. 

The only program changes required would be the addition of the 
appropriate blade property's computer name to the NAMELIST /IN/. 

Since the NAMELIST is read after the blade properties are defined by 
the presently implemented read statement, the quantities read from 
the NAMELIST would replace any previously defined quantities. 
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Program Limitations 


Program limitations fall into two general categories, those 
related to computer compatibility and those related to limitations 
of the model and methods of solution. Some of these program limi- 
tations have been dealt with in other sections, and are repeated 
here for emphasis and convenience to the user. 

Computer compatibility.- One of the most serious limitations of 
the program as listed in APPENDIX C is that it is not functional on 
computers which use less than 14 digits in floating point arithmetic 
operations, due to numerical inaccuracies in the determination of 
natural frequencies and mode shapes. That program is compatible 
with CDC 6600 or CDC 6400 computers at NASA/Langley. Minor changes 
would be required to make the program compatible with other computers 
which have at least 14 digit words. 

Program changes which resulted in satisfactory program operation 
on an IBM 360/65, are given in detail in APPENDIX E. 

Model limitations. - There are basic limitations of the model or 
of the computer program as implemented which restrict the types of 
systems which may be analyzed by the program. The lumped parameter 
nature of the model requires the assumption of uniform property 
massless elastic sections and point masses and inertias; stepwise 
blade twist, and stepwise elastic axis and center of mass location 
variations, with the elastic axes of all elements assumed to be 
parallel to a radial coordinate. Chordwise offsets in the elastic 
axis are assumed to be rigid, but may be large. Shear deformations 
are neglected, so short, thick beams or other systems for which 
shear deflections may be important as compared to bending deflections 
would not be accurately represented by this system. Oscillatory 
axial force variations and displacements are not included as part of 
the model, nor rotary inertia about the z-axis. The blade elements 
may either have an axial tension load as determined by the sum of 
the centrifugal forces on the masses outboard of the blade element, 
or may have zero axial load. They may not be in compression, and no 
spanwise distribution of axial tension load is allowed, other than 
that which is a stepwise (uniform within each element) approximation 
to a centrifugal loading distribution. Only those hub and tip 
boundary conditions which are discussed in the section on PROGRAM 
IMPLEMENTATION are presently available; others would require program 
changes. The hub and tip boundary conditions are assumed to be 
uniform azimuthally. A perfectly rigid hub is assumed with.no azi- 
muthal variations in control system stiffness or bending stiffnesses, 
and no effects of shaft flexibilities or helicopter engine or fuse- 
lage mass are included. 
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Numerical Accuracy Problems 


Numerical accuracy problems are inherent in most eigenvalue 
problems of this type, but have been largely overcome by the modi- 
fied transfer-matrix approach which is implemented in this program. 
There are two general areas where numerical inaccuracies may be a 
problem. One is the use of the frequency search technique to obtain 
natural frequencies and the other is the definition of acceptable 
mode shapes with a convergence criterion based on the smallness of 
the change in the estimated natural frequencies for two successive 

A 4- n 4 — i r\nc* 
uur u u-LWiio • 


Effects of numerical inaccuracies in the frequency search 
options.- The types of numerical accuracy' problems which are associ- 
ated with the frequency search methods have been reduced by the auto- 
mated search method, but are still real. The use in the iteration 
process of one of the remainders rather than the determinant or 
auxiliary function may result in convergence to a frequency other 
than that intended, regardless of the type of search method used. 

That is, duplicate frequencies may be obtained and frequency skip- 
ping may occur. This is less likely with the automated search option 
than with either of the other methods. A discussion of other types 
of problems which may be encountered with the frequency search 
options is given in the THEORETICAL FORMULATION section. 

The frequency iterations are now stopped by either a limit on 
the number of iterations or convergence of the frequency to a natu- 
ral frequency value. If the iteration limit is reached, the program 
will print out a mode shape corresponding to that frequency. Usually, 
this mode shape does not satisfy the tip boundary conditions. A 
satisfactory mode shape may be obtained in most cases by a subsequent 
run using the initial estimate option with an initial estimate as 
obtained from the output of the run which failed to converge. When 
the automated search method is limited by the number of iterations, 

the natural frequency counter, NFREQ, is set equal to NUMF , and the 
program will go to the next model. This is necessary to avoid in- 
correct definition of the auxiliary function and subsequent erratic 
behavior of the automated frequency search. Thus, a higher value for 
NMAX is suggested for use with the automated frequency search option 
than with the other options. 

Effects of numerical inaccuracies on mode shapes .- The use of a 
convergence criterion which depends on the smallness of the change 
in successive estimated natural frequency values may not always re- 
sult in satisfactory mode shapes, as evaluated by the relative 
orders of magnitude of tip force type quantities as compared to 
these quantities at the next inboard section. However, use of 

NEXPO = 12 has resulted in satisfactory mode shapes for all cases 

which have been run with this program, and use of NEXPO = 10 is 

usually satisfactory. It is expected that use of NEXPO > 13 is futile 

for computers which use 14 digit floating-point words. 
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DISCUSSION OF TEST CASES 


Test cases were run at NASA Langley which tested all computer 
program subroutines and which included all program control options 
available for frequency search, mode coupling, boundary conditions, 
and punched output (including suppression of punched output of tor- 
sional mode data) . The automated frequency search method was used 
with helicopter blade models which had been made with the previous 
program to insure consistency of results, and for polynomials with 
multiple and closely spaced roots to provide a test of the automated 
search method for possible similarly severe cases which may be en- 
countered in future use. All options perform as expected, and the 
automated search technique handles multiple and closely spaced roots 
very well. Details of some of the test case results are given 
subsequently. A listing of the data deck used to generate the test 
cases is presented in APPENDIX D, together with selected printed 
and punched output. The complete printed output has been provided 
to NASA Langley, and may be duplicated by running the program with 
the test case data deck. This run was made during March 23, 1972, 
and required approximately 168 CPU seconds (including 13.4 CPU 
seconds to load and compile) and 65 PPU seconds for loading, compil- 
ation, and execution. The storage required as given by the FWA 
LOADER and LWA LOAD labels was 041333 and 032750 respectively; 

276 O/S CALLS were made during the run; and 5120 lines of printed 
output were generated. 
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Test Cases For Helicopter Rotors 

Helicopter rotor test cases included runs for an approximately 
uniform rotor with various boundary conditions, a fully articulated 
(hinged-hinged) rotor with various types of mode coupling, and one 
with two nearly coincident frequencies. All frequency options were 
tested and operated as expected. Most of the cases chosen for dis- 
cussion used the automated search method, but one frequency stepping 
case is also presented. 

Various boundary conditions .- Test cases for helicopter rotors 
which had been run with the previous program included that labeled 
in APPENDIX D as Test Case l.a. Approximately Uniform Clamped- 
Clamped Rotor. Three frequencies were determined, and agreed with 
those determined by previous runs (not using the automated frequency) 
search method). Those three frequencies were 1.0452ft, 1,5110ft, and 
2.7756ft. No printed or punched mode shapes corresponding to these 
frequencies are included in APPENDIX D. The first frequency and 
mode shape for this same blade but with pinned flatwise-clamped 
edgewise (IBC=3) and clamped flatwise-pinned edgewise (IBC=4) 
boundary conditions are presented in APPENDIX D as Test Case l.b 
and l.c, respectively. 

Mode coupling variations .- All types of mode coupling were ex- 
ercised for a fully articulated rotor. Of these, two are presented 
in APPENDIX D as Test Cases 2. a and., 2.b. The rotor parameters and 
program control parameters are given with Case 2. a, which is for the 
first frequency for uncoupled torsional motion (NTYPE=1) . Only NTYPE 
is different for Case 2.b, which gives the third frequency for par- 
tially coupled torsional-f latwise motion (NTYPE=4) , and is the par- 
tially coupled equivalent of the first uncoupled torsional mode. 

Closely spaced frequencies .- A model of a fully articulated 
rotor with very little pitch- flap coupling was modified to give two 
fully coupled (predominantly torsional-f latwise) modes with nearly 
coincident frequencies. Test Cases 3. a and 3.b present the blade 
model parameters, program control parameters, frequency search data, 
and mode shapes for the two coincident frequencies. Note that the 
option to eliminate previously determined natural frequencies was 
also used during this run. 

Frequency stepping case .- The frequency search data presented 
in APPENDIX D as Test Case 4 indicates the type of problems which 
may be encountered if the frequency stepping option is used. The 
blade data for this case is the same as that in Test Case 3. It 
should be noted that the choice of N00R=1 was used, which results in 
a torsional-motion dominated search. Use of other remainders (2,3, 
or 4) may have resulted in a search which did not miss the frequen- 
cies near 53.95 rad/sec and 70.67 rad/sec, as occurred with the use 
of NOOR=l . 
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Polynomial Case 


For a model with NSEC=2, the program automatically assumes that 
the squares of NUMFeNPDF roots of a polynomial are read from the 
OMINT array. A polynomial with roots at the square roots of the numbers 
5, 6, 6, 10, 10.01, 23.47, 33.8, 999, 1001, and 15000 was used to 
test the automated frequency search option. In this operation, the 
numerical value of the polynomial was determined by the program, 
and that value was treated as a determinant for an uncoupled tor- 
sional mode type. The table below gives the squares of the actual 
and the computed roots as found using the automated search method. 

As may be seen, the method successfully obtained the multiple root 

at / 6 , and also distinguished between the closely spaced roots at 

/lO and /10 . 01; and at /999 and /10Q1. 


Values of Squares of Roots 


actual 

computed 

5 

4.9999999997 

6 

5.99999999998 

6 

5.99999999997 

10 

9.999999997 

10.01 

10.00999996 

23.47 

23.46999997 

33.8 

33.7999998 

999 

998.999999997 

1001 

1000.999998 

15000 

14999.547 
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k 


part of 
A , AS AV 


A 


IXMBOLS 


:h are used in this document are given 
ilent, if any, and definitions. 


O T TTYlK/^ 1 
mis-/ w _i_ 


i/EI , length flexibility parameter, 
1/foot-pound 

k^, as determined from "excluded 

equation" of equation (10) ; for 
example, a^ as given by equation (12) 

k^, as determined from simultaneous 
solution of equation (11) 


cosh y , see APPENDIX B 

(cosh y- 1)/ 2 / see APPENDIX B 

rigid offset matrix, see equation (2) 
and APPENDIX B 

de terminant 

elastic field matrix, see equation 
(2) and APPENDIX B 

torsional, edgewise and flatwise 
stiffness, respectively, pound-feet^ 

i^* 1 auxiliary function 
frequency step factor 

offset of elastic axis ahead of radial 
coordinate, in plane perpendicular to 
shaft, feet 

iteration counter during frequency 
search 

rotary inertia about x and y axis 

2 

respectively, pounds feet seconds 

mode shape correction coefficient 
matrix, see equation (10) 

torsional control system stiffness, 
foot-pounds/radian 
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Symbol 

Algebraic 


Name (s) 
Computer 


Symbol Definition 


l 

l 

a 


c 

z 


M 

Y 



m or m 


N 


[P] 



[Q] 

q 

R 

r 



v 

w 

X 

root 
x ,y , z 

(z) 

Y 


ELNTH 

ELC 

ELZ 

Z (9) ,Z (5) 


EMAS or AVMAS 
NR 

PRD , SAV 
PTX,PTZ 

SAV 

NFREQ 

ERR 
R or X 
S 
S 
T 

z (2) 

-Z (6) f Z (10) 

Z(3) 

Z (7) 

XROOT 

x,-,- 

z 

GAMMA 


EPS 


section length, feet 

control link length, feet 

offset of elastic axis in chordwise 
direction, ft 

z axes, respectively, foot pounds 

point mass and inertia matrix, see 
APPENDIX B 

element mass or lumped parameter 

2 

point mass, pounds-second /foot 

number of elements in the state 
vector, (z) 

product matrix, see equation (6) 

centrifugal forces is the x and z 
directions, respectively, pounds 

determinant matrix, see equation (7) 

number of frequencies found, see 
figure 8 

remainder 

radial coordinate, feet 

sinh y/Y/ see APPENDIX B 

(sinh y~y)y 3 / see APPENDIX B 

transfer matrix, see equation (2) and 
APPENDIX B 

torque, foot pounds 

shears in the y. and z directions, 
respectively, pounds 

flatwise lineal deflection, feet 

edgewise lineal deflection, feet 

length from center of rotation to 
inboard end of first section, feet 

axial, flatwise, and edgewise local 
blade coordinates, see figure 1 

state vector, see equation (1) 

MP t /EI) 1/2 , see APPENDIX B 
x 

mass eccentricity, see APPENDIX B 
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Symbol 

Algebraic 

Name (s ) 
Computer 

Symbol Definition 

0,4., ijj 

z (4) ,z (1) ,z (8) 

angular elastic deformations about 
z, x, and y axes, respectively, 
radians 

Ik] 

ERR 

state vector correction factor matrix 

IX] 

n. 

B 

estimated hub state vector quantities, 
see equation (8) 

[A] 

part of 
A , AS AV 

estimated tip state vector, quanti- 
ties, see equation (9) 

In] 

A or ASAV 

modified matrix transfer matrix 

[*] 


twist transfer matrix, see equation 
(2) and APPENDIX B 

$ 

CAPHI 

angle between local blade element 
z-axis and plane perpendicular to 
shaft, degrees for input, radians 
internally 

Q 

CPOMG 

rotor speed, radians/second 

0 ) 

OMINT 

frequency estimate, radians/second 

w n 

OMEGA 

natural frequency of lumped parmater 
model, radians/second 

Superscripts 



i 


iteration number 

Subscripts 



a 


actual value 

e 


effective value 

i,k 


matrix element row and column position 
indices, respectively 

j 


blade element number 

max 


maximum value of parameter 

N 


tip end quantity 

0 


hub end quantity 

x,y,z 


local blade coordinates parallel to 
the radial coordinate and in the flat- 
wise and edgewise directions respect- 
ively, see figure 1 

Mathemtical 

Notations 


A 


an incremental change or amount 

0 


d_ 

dt 
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APPENDIX B 


TRANSFER MATRICES 

The transfer matrices which are used in defining the element 
transfer matrix, [T] , of equation (2), are given here for conveni- 
ence. The computer program implementation was done on the basis of 
the analytically defined products, i.e. the elements of [T] are 
programmed. For convenience, these matrices are partitioned into 
submatrices corresponding to uncoupled mode state vectors. 

Torsional twist matrix, [$], for an incremental twist angle $ , 

0__ __ _ l_ _ 0 __ _ _____ __ 

cos $ 0 0 0 1 -sin$ 0 0 0 

0 cos$ 0 0,0 -sin$ 0 0 

0 0 cos$ 00 0 -sin$ 0 

0 0 0 cos $ 10 0 0 -sin$ 

sin<l> 0 0 0 | cos $ 0 0 0 

0 sin$ 000 cos$ 0 0 

0 0 sin$ 0 I 0 0 cos$ 0 

0 0 0 sin$ | 0 0 0 cos4> 

Massless elastic field with tensile force matrix, [E] : 
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where a 


= */EI, 



cosh y 

2 

Ccosh y - 1}/y 
sinh y/y 


Csinh y ~ y)/y" 


P T /El 

X 


>1/2 


p m = « Z rn.r. 
x i=j 

where j is the index of the current element being consi- 
dered and m. is the lumped mass at a radial position of 

r i* 

The subscripts x, y, and z refer to the El, S^, S C^, and C 2 

quantities within the appropriate brackets. That is the (2,2) ele- 
ment is and the (4,5) element is 1(1/ EI Z ) (sinh Y z /y z ) • 


Point mass and inertia (in a centrifugal field) matrix, {M] : 


1 

0 | 

0 

0 

0 

0 1 

0 

0 

0 

0 

J*2,l 

1 

| M 2,3 

0 

0 

! 

M 2,7 

M 

2,8 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 ^ 

0 

0 

0 

0 

M 5,l 

0 ^ 

0 

0 

1 

0 

0 

M 5 , 8 

0 

0 


0 

Vs_ 

0 

0 

1 

J 

'l 

^6,8_ 

0 

0 

0 

0 

0 

0 

0 

0 1 

0 

0 

0 

0 

0 

0 

0 

0 

0 ! 

' 0 

1 

0 

0 

M 9 , 1 

0 

M 9,3 

0 

0 

0 

1 M 9,7 

M 9 , 8 

1 

0 

M 10,l 

0 

M 10,3 

0 

0 

0 

1 M 10 , 7 

M 10 , 8 

0 

1 
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where M 


ft ) + hemft C 


= I x a + lyB CC^, - SJ1 + £ m Ca - 
M 2 3 = -em (a 2 - ft 2 S 2 ) 

M 2 , 7 = “ emfi S y c <i , 

M 2,8* = " fia (I x + r y + e2m) S * 

2 

M . = remft 

b / X 

M g g = -ft 2 [CI X + e 2 m) + hem] 

M 6 ^ = em (a 2 - ft 2 S 2 ) 

M g ^ 3 = -in (a 2 - ft 2 sj) 

M 6,7 “ -" a2 S t S 

Mg^g* = -fto (2em) 

M 9,l* = {I x + X y + e2m) S ¥ 

Mg g* = -fto (2em) 

Mg y* = fta (2em) 

Mg 8 = “I^S 2 + I 0 2 + e 2 m ( a 2 - ft 2 S*) + h e m^ 2 C,,, 



M 10 , 3 

2 

— S C 

M 10 , 7 

-m (c 2 - a 2 C 2 } 

II 

* 

CO 

K. 

0 
1 — 1 
s 

Ha (2em) C 

5 

■> 

11 

1 


and all other M. . = 0 

1 r 3 

Here m is the lumped mass, and is the average of adjacent sec- 
tion masses, and e is the corresponding averaged eccentricity. 


a indicates -57- 

dt 

= sin f 
Cy = cos V 

and ¥ is the angle between the load blade element, and the plane 
perpendicular to the rotor shaft. 

Terms with a * superscript are gyroscopic terms, and are not 
used in the program as presently implemented. (Their use would 
result in complex mode shapes, which are not orthogonal in the 
usual sense.) In any blade loads and response program using purely 
real modes as generalized coordinates, these gyroscopic terms could 
be used as part of the forcing function. 
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APPENDIX C 


PROGRAM LISTING 


Computer program listings are given with each subroutine be- 
ginning on a new page. Job control cards are not included, as these 
vary from one facility to another The program does not use OVERLAY 
or SEGMENT options. Program output format statements give quantities 
labeled in English units. Input may be in any consistent set of 
units, but output would be incorrectly labeled unless: (1) English 
unit input quantities are used or (2) output formats are changed. 

The listing as given here is identical to that used on CDC 6600 
computers at both the NASA Langley facility and at a non-governmental 
facility. The program should not be expected to operate properly on 
computers which use fewer than 14 digits for floating point numerical 
operations . 

The main program and subroutines are listed in the approximate 
order of their first use in the program, with page numbers noted for 
convenience, as follows. 

BLADE 6 6 

SNTFR 7 5 

BNDRY1 76 
BNDRYX 81 
TEMAT 8 4 
MATINV 93 
FINDM 95 
INTERP 98 
SHAPE 99 

Program structure is indicated by the following chart. 



65 




o o 


PROGRAM BLADE! I NPUT* 1001, OUTPUT* 1001 * TAP E5* INPUT ,TAPE6=0UTPUT» 
IPUNCH*1001 *TAPE4* PUNCH, TAPE 7« 1001 • TAPE 1 , TAPE99 ) 

C MAIN PROGRAM FOR MATRIX METHOD SOLUTION OF DYNAMIC BLADE PROBLEM 
DIMENSION ELI30), EIXI30)«E1VI 301 ,EI 2 (30) • XINRI 30 1 , YINRI 30 ) 
DIMENSION EMASI 30),PH1(30),EPSI30),ELZ(31),CCPS(30) , CAPS I 130 ) 
DIMENSION PROI 10, 51 ,S AVI 10, SI ,T( 10, 10) ,ASAV( 10,51 ,AVMASI3ll 
DIMENSION ROLOI4) ,RNEM(4) ,ERR(10),A( 10,51 ,B( 10) , Cl 101 
DIMENSION NANDI 5 1 ,GM!NT(15) 

DIMENSION TITLE 119) 

DIMENSION DATEI2) ,PHI XI 30) , EPSXI 30) 

DIMENSION C101I 10,1) ,ERR101 (10,1 ),SV1010( 10,10) 

COMMON CPSQ«GMSQ,Xkuui,£EKMG 

COMMON EL,EIX,EIY,EIZ,XINR,VINR,EMAS,PHI «EPS ,ELZ,CAPSI ,CCPS, EPSN 
COMMON CP0MG,0MEG,S1ZER,THETC,ELC, AVMAS 
COMMON NGYRO.NX ,NY,NZ , NSEC, NS PI 
COMMON NC ,NC1 ,NC2 ,NC3 «NR,NR1, NR2 «NR3 «NR4 
COMMON NIBC,B 
COMMON ERR,TITLE,NFREQ 

EQUIVALENCE IC101 1 1 ) »CI 1 ) ) • IERR101I 1 ) ,ERRI 1) ) , I SV1010U ) , SAVU ) ) 
NAMELIST/ IN/NTVPE,NMAX,NEXPO, NUMF,NOOR,NPCH,NFINO,NPDF,CPOMG, 
1SIZER,THETC,ELC,XR00T ,PRCNT ,OMINT ,TMSL,OMLIMT 
SEE PAGE 3 FOR GEOMETRY 

CALL OAYT IMIDATE) 

N0UT«6 
C INPUT 

IBCP-O 

420 READl 5,401 ) IBC, I TAB, INSEC, TI TLE 
IFIINSEC.GT.30) GO TO 1200 

C REPLACE THE NEXT CARD BY A TEST ON END OF DATA OR REMOVE 

430 IFIITAB.EQ.O) GO TO 440 
NSECMNSEC 

IFINSEC.NE.O) GO TO 433 
PRINT 431 

431 FORMAT I33H NSEC-O, CANNOT READ IN NEW TABLE ) 

STOP 

433 DO 435 I*1,NSEC 

REAOI 5,402) J*EL( J) ,E IXI J1 , El Y< J) ,E I Zl J) , XINRI J) , YINRI J ) 
l.EMASI J),PH1( J),EPS(J),ELZI J) 

PHIXI J)«PHII J) 

EPSXI JI-EPSIJI 
435 CONTINUE 

EPSN«EPSINSEC) 

440 REAOI 5,IN)- 

IFINFIND.GE.O) GO TO 338 
NSTEP-1 
NTAPE»-i 
GO TO 339 
338 NSTEP«0 
NTAPE*0 
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339 CONTINUE 
NTRL*NPDF 

I F( NF IND.GT *0 1 NFLG-0 

C DEFAULT FOR NTVPE OUT OF ALLOWABLE RANGE 

4401 IFINTVPE.GT.O.ANO.NTVPE.LT.?) GO TO 4402 
CALL SNTFR I NTYPE*NR1*NR2*NX*NV*NZ* NBNO) 
GO TO 430 

4402 GO TO (441*442 *443*444*445*446 ) t NTYPE 

441 NR1*1 
NR2*2 
NX* 0 
NY* 1 
NZ* 1 
NBNOI 1 1*2 
GO TO 450 

442 NR1*3 
NR2*6 
NX *1 
NY *1 
NZ *0 
NBNOI 1)*5 
NBNOI 21*6 
GO TO 450 

443 NR1*7 
NR2*10 
NX *1 
NY *0 
NZ *1 
NBNOI 11*9 
NBNOI 21*10 
GO TO 450 

444 NR1*1 
NR2*6 
NX* 0 
NY* 1 
NZ* 0 
NBNOI 11*2 
NBNOI 21*5 
NBNOI 3)«6 
GO TO 450 

445 NR 1*3 
NR2*10 
NX* 1 
NY* 0 
NZ* 0 
NBNOI 11*5 
NBNOI 2)*& 

NBNOI 31*9 
NBNDI4)*10 
GO TO 450 
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446 NR1*1 
NR2*10 
NX* 0 
NY* 0 
NZ* 0 
NBNOI 11*2 
NBNOI 3)*6 
NBNOI 21*5 
NBNOI 41*9 
NBNOI 51*10 
450 0NEG*0H!NTI1I 
SIZERX*$iZER 
THETCX*TMETC 
NbYR0*0 
NCl*t NR1+1 1/2 
NC2*NR2/2 
NC*NC2*NC1*1 
NC3*NC2-NC1 
NC4-NC-1 
NC5-NC-2 
NR»NR2— NR141 
NR3*NR-1 
NR4*NR-2 

1FIN00R I 5123,5123*5122 

5122 IFIN00R-NC4) 5124,5 124,5123 

5123 N00R*NC4 

5124 CONTINUE 
WRITEINOUT,955) 

PRINTOUT OF INPUT 
WRITE INOUT, 950) 

WRITE 16,8981 TITLE, 0ATEI1 » 

C SPECIAL MESSAbES IF INFINITE STIFFNESSES OCCUR IN ANY OIRECTION 
481 FORMATI 16151 

C LABEL THIS WRITE IF TIME PERMITS 
75 IF I NX) 80,80,77 
77 WRITE INOUT ,956) 

80 IF INYI 85,85,83 
83 WRITE INOUT ,957 ) 

85 IF INZ) 90,90,87 
87 WRITE INOUT, 958) 

90 WRITE IN0UT,96l)NSEC,CP0M6,SIZER,THETC,ELC,XR00T 
WRITE INOUT ,962 ) 

WRITE IN0UT,965) II ,ELII ) ,EMASI I),EIXtl),EIY( I) ,EIZ( I )« XINR( I ) , 
1 YI NR II), PH II I),EPSU) fELZII) ,1*1, NSEC) 

192 CONTINUE 

WRITE INOUT, 955) 

WRITE INOUT, 403) NTYPE.IBC .NFIND, I TAB , TMSL ,NUMF 
1 , NPCH,NP0F,NMAX , NOOR, NEXPO , CMLIMT , PRC NT 
IFINFIND.EQ.1.AND.NPDF.GT .0 ) WRITE! 6,405 ) I OM INTI I ) , I = 1,NPDF) 
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IFINF IND*EQ*O.AND.NUMF.GT .0 I WRITE! 6*4041 (OMINTI I ) » I=1,NUMF I 
1FINFIND«LE*0 .AND.NUMF.EQ *0 I WRI TE( 6*4031 1 
WRITE (6*9811 

IFINC4.GT.0I WRITE (6*982) (N«N»1*NC4) 

IFINC4.EQ.0) WRITE (6*982) 

DEFINE CONSTANTS CONVERT INPUT DEGREES TO RADIANS 
NSP1-NSECM 
NFREO-l 

IFINF IND.LE.O) GO TO 36 
NFREQ-NPDF 

iFi NFREQ.GT.O ) GO TO 36 
OMEG*l. 

GO TO 37 

36 OMEG-OMINTI1) 

37 PI»3* 1413926 
CNVRT-PI/180. 

DET-0. 

IF (NFIND.GT.O) QMSQ»1. 

NTIME-O 
NFLAG-0 
RN-0. 

IF (NSEC.EG.2) GO TO 277 

SIZER«SIZERX*CNVRT 
THETC-TMETCX*CNVRT 
C PSO*CPONG*C PONG 
OHSQ«OMEG*OMEG 
CAPSI (11-SIZER 
CCPSI D-COS1SIZER) 

39 00 40 I-2*NSEC 
INl-I-i 

EPSII )-t EPSXI 1 )*EMAS( II*EPSX( IM1 )«EMAS(IM1) )/t EMASt I )+EMAS( I HI ) ) 
PHI ( I )-PHIX( I )*CNVRT 
CAPSI (I)-CAPSK I-ll+PHI (I » 

40 CCPSI I l-COSICAPSI ( I ) ) 

41 IF (PRCNT) 60*50*60 
50 PRCNT-10. 

60 FACTR-1.+PRCNT/100. 

C INITIALIZATION 

65 DET-O* 

IFINF INO.GT.O) OMSQ-1. 

NTINE-0 

NFLAG-0 

RN«0. 

DO 92 N«l *4 
B(N)*1. 

CINI-O. 

92 RNEWt N) s 0. 
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Bt5l*l. 

Ct5l*0. 

DEFINE 8QUN0ARV MATRIX 
100 1*1 

IF fNSEC.EQ.2l GO TO 277 

IFfIBC.EQ.il GO TO 471 
CALL BNORVXf PRO. A. I • IBC) 

GO TO 475 

471 CALL BNDRYlf PRD.A.II 
475 NIBC*I 

DO 180 I*NIBC.NSP1 
00 101 , 1 * 1.10 
DO 101 JX-1.5 
ASAVf J.JXI-O. 

101 SAVf J « JXI*0. 

00 281 J*1.10 
DO 281 JX*1.10 
281 T f J«JXI*0. 

CALL TENATfl.T) 

MODIFY BNDRY SUBROUTINE FOR ANY CHANGE IN BOUNDARY CONDITION AT 
CENTER ENO 

PREMULTIPLY BY NEW T MATRIX 
DO 140 M*NR1.NR2 
DO 140 N*NC1.NC2 
DO 140 NINT*NR1«NR2 

ASAVf M.NI*TfM.NINTl*AtNINT,N> ♦ASAVf M.N) 

140 SAVfM,NI*?fM,NINTI*PRDfNINT,N)*$AVfM,NI 
DO 180 M«NR1.NR2 
DO 180 N*NC1.NC2 
At M.N I* ASA VI M.N I 
PRDtM.NI*SAVf M.NI 
180 CONTINUE 
200 CONTINUE 


COMPUTE SHAPES FOR ITERATION ON FREQUENCIES 

EXTRACT SUBMATRIX FOR BOUNOARY CONDITIONS AT THE TIP 

BOUNDARY CONDITIONS AT TIP ARE CONTAINED IN EQUATIONS 2*5,6.9.10 

GIVING 5X5 SYSTEM TO FINO EIGEN VALUES OF 

IFtNR2>2 1 277,277, 249 
249 DO 250 M*1.NC 
J*NBNOf Ml 
JKL=NCl-l 
DO 250 N* 1 » NC 
JKL=JKL+l 
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ASAV<M,N)*AIJ,JKL) 

250 SAVIM,N)*PRDIJ,JKL) 

NXl»l 

00 276 J*l,NC4 
1*0 

00 254 M*1,NC 
IF IM-JI 252*254,252 
252 1*141 

ERRII)*-ASAVIM,1) 

00 253 N*1,NC4 
T(I«NI»A$AVIM,N4l) 

26 3 CONTINUE 
254 CONTINUE 

CALL MATINVI T,NC4, ERR 101 ,NX1, OT, 101 1 
IFIABSIDTI-. IE- 201256, 257,257 

256 URITEINQUT,980) 

257 ROLDIJ)*RNEUtJI 

I FI ABSI ASAVIJ.J41 » I-. IE-201 258,258, 259 

258 RNEUI J )*Q. 

82*0. 

GO TO 276 

259 B2«-ASAVIJ,1) 

00 270 N*i*NC4 

1 FI J-N) 265,270, 265 

265 B2*B2-ASAVI«J,N4l ) *ERR INI 
270 CONTINUE 

B2*B2/ASAVU,J41) 

RNEUI J)*ERRt JI-B2 

275 CONTINUE 

272 IFI J-N00RI276 ,273,276 

273 CUI*BIJ|4.5«IERRI J)4B2) 

00 274 LL*1«NC4 

274 aiLL)*B<LL»4ERR(LL) 

BIN00R)*C( II 

276 CONTINUE 

ENTER MATINV TO 

COMPUTE DETERMINANT VALUE OF 5X5 SYSTEM 
OTOLO*OET 

CALL MATINVI SV1010,NC,C101 ,NX1,DET, 101) 

OMEGA*SQRT I OMSO ) 

IFINF INO.GT.O ) GO TO 721 

URITEINOUT,972) NTIME , OMEGA ,DET, ( RNEUI J) ,J*1,NC4) 

GO TO 3278 

721 CALL FINOMIOMLST,OMSQ,OET,OTOLO,RNEW,ROLO,NEXPO,NC4,NFREO,OMINT,NT 
IIME «NFIND,NTRL*NFLG«NPDF) 

GO TO 1500,300) ,NTRL 

277 0T0LD*0ET 
DET*PRDI 2, 1 ) 
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IF (NSEC. NE. 2 1 GO TO 2772 
DET*l.O 

00 2771 I»1.NUMF 

2771 DET*(OMINT(I)-OMSQ)*OET 

2772 CONTINUE 
C 

IF(NF 1N0.GT.0) GO TO 721 
N00R*1 

RNEM(NOOR)*DET 
ROLDI NOOR) »DT OLD 
3278 IF INSTEP) 279.279.278 
C 

278 IF (NFLAG) 2781.2781.285 

C CHECK FOR CHANGE IN SIGN OF DETERMINANT FROM PREVIOUS VALUE 

2781 IF(DTOLD*DET) 2782.280 .280 

2782 NFLAG»1 
GO TO 285 

279 IF (NTIME) 280.280.285 

280 OMLST-ONSQ 
OMSO»OMSO*FACTR 
GO TO 290 

285 I F(NR2'*’2 1286.286. 287 

286 RNEM(NOOR)«DET 
ROLDI NOOR) *DTOLO 

C INTERPOLATE LINEARLY FOR NEW VALUE OF OMSQ 

287 YAsRNE Ml NOOR ) 

YB*ROLD( NOOR) 

CALL INTERPIOMLST .OMSQ. YA.YB) 

290 NTIME*NTIME+1 

C CHECK FOR REQUIREO CONVERGENCE ON VALUE OF OMEGA SQUAREO 
IFUBSI0NSQ-0MLST)/0NLST-.1**NEXP0) 500.500.300 
300 IF (NTIME.LT.NMAX ) GO TO 100 
MRITE (6.984) 

IFINF INO.GE.l ) NFREQ-NUMF 

C COMPUTE DISPLACEMENTS ANO FORCES ALONG THE BLADE 
500 MRITE (NOUT. 955) 

MRITE (6.898) TITLE. DATEd ) 

NTIME*0 

IF (NSEC.EQ.2 ) GO TO 600 

CALL SHAPE ( PRO. NT APE • NPCH. I BC . TNSL ) 

IF (NUMF) 800.800.600 
600 IF INFREQ-NUMF) 6C5.800.800 
605 IF ( OMINT I NFREQ ) • GT .OML IMT ) GO TO 800 
NFRE0*NFRE0^1 
MRITE (6.983) 

IF(NC4.LE.O) MRITE (NOUT. 982) 

IF(NC4.GT.O) MRITE (NOUT. 982) (J»J*1»NC4) 

607 IF (NFINO.GT.O) GO TO 721 
IF(NFINO.GT.O) GO TO 721 
OMSQ*OMSQ*FACTR 
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IF !NTAPE) 620 *610*610 

C ALLOWS FOR NEW INITIAL GUESSES FOR EACH FREQUENCY UNLESS NTAPE=-1 
610 OMEG*OMINT!NFREQ) 

OMSQ»OMEG*GMEG 
NSTEP*0 
GO TO 65 
620 NSTEP*1 
GO TO 65 
800 CONTINUE 
1130 CONTINUE 
GO TO 420 
1200 STOP 
C 

401 F0RMATI2I 1*I2«19A4) 

402 FORMAT! 1 2* F6.5.9E8.0) 

898 FORMAT!22X*19A4*l0X*A10*/ I 

950 F0RMAT!50X,23HBLADE FREQUENCY PR0GRAM//38X*46HBLADE MODEL AND OPER 
IATING CONDITION PARAMETERS//) 

403 F0RMAT!48X,26HPR0GRAM CONTROL PARAMETERS//10X, 22HM0DE COUPLING* NT 
1VPE * • I 3 • 10X * 28HHUB BCUNDRY CONDITION* IBC », 13* 15X* 25HFREQUENCY S 
2EARCH* NFIND », I3/5X, 27HBLADE PROPERTY DECK, ITAB =, I3,ilX,27HT0RS 
3ICN MOOE SWITCH* TMSL «,F5. I* 12X,26HN0. OF FREQUENCIES* NUMF =, 
4I3*/10X*22HPUNCHED OUTPUT* NPCH *, 13 ,6X,32HN0. OF KNOWN FREQUENC IE 
5S* NPOF »• I3*15X*25HN0. OF ITERATIONS* NMAX * * 1 3*/9X*23HREMAINDER 
6INDEX* NOOR »*I3, 12X.26HC0NVERGENCE LIMIT* NEXPO »,I3,/7X, 
725HFREQUENCV LIMIT. OMLIMT », E11.3, 5X,25HSTEP SIZE FACTOR, PRCNT = 
8.F11.6) 

404 FORMAT (4X*28HFREQUENCY ESTIMATES, OMINT =,2X, ( 7F12.6) ) 

405 FORMAT! 3X.47HPREVI0USLY DETERMINED FREQUENCIES, OMINT(I) * , ( 7F1 

12 . 6)1 

955 FORMAT ! 1H1I 

956 FORMAT !43X*37HINFINITE STIFFNESS IN THE X DIRECTION) 

957 FORMAT !43X*37HINFINITE STIFFNESS IN THE Y DIRECTION) 

958 FORMAT ! 43X.37HINFINI TE STIFFNESS IN THE Z DIRECTION) 

961 FORMAT! //43X.20HNUMBER OF SECTIONS »,I3,/35X, 

A 28HR0T AT IONAL SPEED CAP 

IOMEGA «F9.4,12H RADIANS/SEC/40X,23HC0NTR0L ANGLE SI ZERO =F9.4, 

2 8H DEGREES/40X* 23HCONTROL ANGLE THETA C -F9.4,8H DEGREES/ 

3 47X « 16HC0NTR0L OFFSET *,F9.4,7H FEET / 23X .40HRAOIAL DISTANCE T 
40 FIRST BLAOE ELEMENT *,F9.4//) 

962 FORMAT I3H I *4X .6HLENGTH, 7X ,4HMASS .9X.3HEI X ,9X*3HE I Y, 9X ,3HEI Z , 

1 8X*2HIX*10X*2HIY*10X* 3HPHI ,7X, 7HEPSI L0N,5X,6H0FFSET /7X , 

2 6H FEET *6X,7HLB-SEC2 ,5X,6HLB-FT2,2< 6X.6HLB-FT2) , IX, 

3 2! 4X, 8HLB-SEC2- ) » 5X ,7HDEGREE S,6X *4HFEET ,8X,4HFEET/19X» 

4 7HPER-FT ,34X ,2 C8X, 4HFEET ) , 2 ! / ) ) 

965 FORMAT II3.10E12.4) 

970 FORMAT ( / *(5E20.12)) 

972 FORMAT! 4X , I3.3X, E20. 12 ,E 15 .5 , I5X.4E 15.5 ) 

980 FORMAT !53X,24H4 X 4 SYSTEM IS SINGULAR) 

981 FORMAT!///, 43X,34HSEARCH FOR LOWEST FREQUENCY BEGINS,/) 
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982 FORMAT! 2X*10HITERAT IONS *5 X*9HFREQUENCY *6X*11HDETERMINANT* 

16X# 9HAUXILIARY *5X* 10HREMAI NDERS* /2X » 6HNUHBER.8X. 9HI RAO/SECI *24X 
2* 0HF UNCTION* 113*3115) 

983 FORMAT! 1H1*47X* 32HSEARCH FOR NEXT FREQUENCY BEGINS*/) 

984 FORMAT!// *23H WARNINGS NTIME » NMAX*//> 

4031 FORMAT! 1X*27HN0 FREQUENCY ESTIMATES USEO) 

7123 FORMAT! I5*5X* 15) 

7124 FORMAT ! 10E8« 1 ) 

ENO 


i 

I 


[ 

I 

t 
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SUBROUTINE SNTFR (NTYPE»NRi *NR2»NX» NY*NZ*NBND) 
DIMENSION NBNDI5I 

DEFAULT FOR NTVPE OUT OF ALLOWABLE RANGE 

GIVES NTVPE*6 CONTROL* I.E.* FULLY COUPLED MODES 
NRl*l 
NR2-10 
NX*0 
NY*0 
NZ»0 

N8N0I1I-2 
NBNDI 21*5 
NBNDI 31*6 
NBNDI 41*9 
NBNDI 51*10 
RETURN 
END 
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SUBROUTINE BN0RY1 (BND.BPL.NPI B) 

SPECIAL SUBROUTINE TO COMPUTE BOUNDARY CONDITION MATRIX 
CASE OP FLAP HINGE.LEAO-LAGHINGE.TORSIONAL SPRING 

DIMENSION SAVXX470) 

DIMENSION BAND! 10*5) 

DIMENSION BND( 10* 5) «BPL 4 10,5).B(10)*SAV( 10*71 

DIMENSION EL(30),EIX430).EIY4 30 ) , El Z 4 30 ) , XINR4 30) * YI NRC 30 ) 

DIMENSION EMAS4 30)«PHI(30),EPS430),ELZ431)«CAPSI430) 

DIMENSION CCPS430I.AVMAS431) 

DIMENSION TITLE (19) • ERR 4 10 I 
COMMON CPSQ.OMSQ, XROOT , ZERMO 

COMMON EL.EIX.EIV.EIZ.XINR. YINR.EMAS .PHI .EPS .ELZ.CAPSI ,CC PS , EPSN 

COMMON CPOMG.OMEG.SIZER.THETC .ELC.AVMAS 

COMMON NGYRO.NX .NY. NZ.NSEC, NS PI 

COMMON NC «NC1 ,NC2 .NC3.NR.NR1, NR2.NR3.NR4 

COMMON NIBC.B 

COMMON ERR. TITLE, NFREQ 

EQUIVALENCE I SAV4 II .SAVXX41)) 

NOUT-6 

GO TO (1.170.201) .NPIB 

1 DO 311 1*1.10 
DO 312 11*1.5 
BPL4 1 .1 1 1*0* 

312 BND4 I * 1 1 1*0* 

DO 313 11*1.7 

313 SAV4 I .II 1*0. 

311 CONTINUE 

L*1 

BND(2.1I*1, 

BND(6,3)*1, 

BND40,4)*1 • 

BND4 10,5 )*1. 

I F( NX ) 2.2,3 

2 BND( 1.1)*EL(1)/EIX(1) 

3 NSM1*NSEC-1 

4 PTX*0. 

X*XROOT 

DO 5 J*1 « NSM1 
X*X*ELIJI 

5 PTX*PTX+.5*4EMAS4 J+D+EMAS4 J) |*X 
PTX-PTX*. 5*EMAS < NSEC) *< X^EL (NSEC ) ) 

PTX*PTX*CPSQ 

IF(NZ) 8,8.45 
8 EI*EIZ(1) 

10 GAMMA«EL41)*SQRT4PTX/EI I 
G2*GAMMA*GAMMA 
IF4G2-.00001ll2.12.il 

11 SOC*( TANH4 GAMMA) ) **2 
S1»SQRT4 SOC/I l.-SOC) l/GAMMA 
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C1*SQRT ( 1. /( 1 ,-SOC) ) 

S2-IS1-1.I/G2 
C2»ICl-i .I/G2 
GO TO 13 

12 Sl*l.+G2/6. 

Cl>l.+G2*.5 

S2»l l.*G2*.05)/6. 

C2-.5+G2/24. 

13 GO TO <20*50)*L 
20 BND<3*2)*EL<1I*S1 

BNDI4,2)*C1 
BNDI 5 «2I*PTX*BND( 3*21 
40 BN0(4*3)*EU1)«EL(1I/EIZ(1) 

BNDI3,3)*EL< 1I*BNDI4,3I*S2 
BND(4*3) *BN0( 4,3I*C2 
BNO( 5 *3) *BN0( 3*2) 

GO TO 48 

43 BND4 3*2) *EL( 1 I 
BNDI4,2)*1. 

BNDI 5,2) S PTX*ELI 1) 

BND45«3)*ELI 1 ) 

48 IFINY) 49*49,55 

49 L-2 
EI-EIVdl 
GO TO 10 

50 BNDI9,5)*-S1/C1 

8ND47»5)*IEL41)**3I/E I*(BND(9 *5) *C2 + S2 ) 
BND49,5)»BN0I9,51*ELI II 
GO TO 58 

55 BN0I9*5I*-ELI1) 

58 IF4ELC-* IE-71 100*100*60 

60 BNDI6«1)*-SINITHETC)/ELC 
BNDI10,1>— COSITHETCI/ELC 

100 ZERM0*BN049*5I 
BNDI9,5>*0. 

00 61 JK*l,10 
00 61 JK1 1*1,5 

61 BANO( JK,JKin»BND(JK,JKII ) 

IFINC-ll 105*105,104 

104 00 101 N 3 1*NC 

101 SAVIN, Nl»l. 

00 102 N*2,NC 

102 SAVIN, 1)*BIN-1) 

105 SAVII,1)>1. 

103 00 111 M«NRi*NR2 
NN«0 

00 110 N 3 NC1«NC2 
NN*NN*1 

110 BPL(N,NCn*BNO(M,N)*SAV(NN,ll+BPHH,NCU 


IF4NC-1) 113* 113*112 

112 NN«NC1*1 

00 111 N»NN*NC2 
111 BPL<M*N)«BN0(M*N) 

113 NPIB*2 

GO TO 400 

170 IF(NC3.EQ.O)GO TO 181 
00 180 N»1,NC3 
180 SAV(N+1,7)*B4N) 

181 X*l. 

BPL410*4)*0. 

BPL<10*5)*0. 

GO TO 203 

201 X«BPL<1*1) 

BPL<10*5)«0.0 

I F4NC3«EQ*0) GO TO 203 
00 202 N*1*NC3 

202 B(N)«SAVf N+l*7) 

203 00 204 II*1*NR 

204 SAVXX(II)*0, 

I»0 

210 1F( NR 1—3 1 211* 212* 213 

211 IF4NR2— 6)216*217*219 

212 IFCNR2-6I224, 224*226 

C EDGEWISE BENOING ONLY NORMALIZE TO W*1.0 

213 V NORM* I BN0(7*4)+BND( 7 , 5)*B4 1) )/X 
SAVI10*l)*B4 D/VNORM 
SAV(9*1)*SAV4 10*1 )*ZERMO 

WRITE ( N0UT*920 )I*4SAV4N*1) ,N*1 , 10) , XROOT 

SAV48»1)*1 • / VNORM 

I>1 

X*XR00T+EL(1) 

C BPL IS RECOMPUTED AS BOUNOARY MATRIX 

214 DO 1214 M*l,10 
1214 B(M)«0. 

00 215 M*7* 10 
00 215 N»4*5 

215 BIM)*BIN)*BAND4 M,N)*SAV(2*N,1 ) 

WRITE 4N0UT. 920)1 ,4 Bl N) ,N*1 ,10) , X 
84 8)*SAV4 8*1 ) 

GO TO 300 

C TORSION ONLY NORMALIZE TO PHI*1. 

216 B(2)*X/BN0( 1*1) 

B 4 1 ) = 0« 

WRITE (NOUT *920)1,641) *B42) , 4SAVXX4II )*II-3*10) , XROOT 

X-XROOT 

GO TO 300 

C FLAP-TORSION NORMALIZE TO V*l. 

217 VNORN*BNDI3,l)*BND43,2)*B< 1 )+BND 4 3, 3 ) *84 2) 

SAV42* i)*X/VNORM 
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SAV(4,l)-B( 1>*SAV<2,1) 

SAV(6,1)«B<2)*$AV<2,1 ) 

WRITE f NOUT * 920 ) I • I SAVI N* 1 1 »N*1 » 10) • XROOT 
00 218 N*l*6 

218 B(N)*SAVIN v l) 

X* XROOT 

GO TO 300 

C TORSION- FLAP-EDGE WISE NORMALIZE TO V*l. 

219 VNQRM*BND1 3,l)+BN0(3»2)*B(l )+BND ( 3, 3 ) *B( 2) +BNDI 3*4 ) *B( 3)+BND(3,5J* 
IBI4I 

SAVI2«1I*X/VN0RM 
SAV(4tl)*Bl ll*SAVf 2«1 ) 

SAVI6«1)*BI2)*SAVI2«1) 

SAVI 10,1 )*B<4I*SAV( 2,1) 

SAVI9,1)*SAVI 10. I )*ZERMO 
WRITE I NOUT 1 920 ) I « I SAVI N« 1 ) ,N*1 , 10) * XROOT 
1*1 

SAVI8*1I»X*BI3I/VN0RM 

220 00 221 M*l 9 10 
BIM)«0. 

DO 221 N*1 *5 

221 B(M|*B(M)*BAN0IM,N)*SAVI2*N,1) 

X*XR00T+ELI1) 

B(8)*SAVI8,1) 

WRITE I NOUT «920 II •(BIN) ,N*1 ,10) ,X 

B PL 110*51*1X1 NR (1 I+AVMASI 2)*EPS(l)*EPS(i) 1*1 Bt 1 )*B( 1 ) +BI 8 )*B (8 > ) 
1-2.4AVHASI2I*EPSI 1 ) *B<3 )*B< 1) ♦AVMAS < 2) *< B( 3) *B( 3 ) +B ( 7 ) *B( 7 ) ) 

GO TO 300 

C FLAPWISE BENDING ONLY NORMALIZE TO V*l. 

224 VNORM*l BNOI 3«2I +BNOI 3*3 )*BI 1) l/X 
SAVI 4* 1 )*1«/VN0RN 
SAV(6«1)*B( ll/VNORM 

WRITE I NOUT *9201 1»(SAV(N»1) *N*1 . 10) • XROOT 
00 225 N»3*6 

225 BIN)*SAVIN«1) 

B(1I*0. 

BI2)*0« 

X* XROOT 
GO TO 300 

C FLAP-EDGE WISE NORMALIZE TO V*l. 

226 VN0RM*BN0I3«2)*BNDI3«3)*BI 1)*BND(3,4)*BI2)+BND<3, 5) *81 3) 

SAVI 4* 1)*X/VN0RM 

SAV(6»1)*BC 1I*SAV(4,1 ) 

SAVI10,1)*B(31*SAV(4,1) 

SAVI 9* 1 )*SAVI 10* ll+ZERMO 
WRITE (NOUT, 920) I* (SAVIN,1)»N= 1,10) .XROOT 
1*1 

X*XROOT*ELI 1 1 
SAV(8,l)*BI2)*SAV(4,l) 

230 00 231 M*1 • 10 


79 



231 B(M)«0. 

00 232 M«3«10 
00 232 N»2*5 

23 2 B(M>*B(M)*BAND(M t N)*SAV<2*Ntl) 
WRITE <NQUT*920tl *(B(N) *N*1*10)»X 
B(8)*SAV( 8*1 ) 

300 NPIB* I 

BPM10«4)*X 
400 CONTINUE 

920 FORMAT (I3*11(E10«3*1XI) 

RETURN 

END 
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SUBROUTINE BNORYX (BNO«BPL *NPI B»KDi) 

OIMENSION BND4 10* 5) *BPL( 10 » 5) .BUOI.SAVC 10*7) 

DIMENSION EL(30)*EIX(30),EIY(30>*EIZ(30) , XINRI 30 I * Y I NRC 30 1 
DIMENSION EMAS( 30 ) • PHI( 30 ) » EPS( 30 ) * ELZ( 31 ) *CAPSI (30 ) 

DIMENSION CCPSC 30 ) • AVMAS( 31 ) 

OIMENSION TITLE(19)*ERR(10) 

COMMON CPSQtOMSQ* XROQT • ZERMO 

COMMON EL*EIX*EIV*EIZ*XINR.YlNRtEMAS.PHI,EPS,ELZ.CAPSI.CCPS,EPSN 

COMMON CPOMG»OMEG*SIZER*THETC *ELC»AVMAS 

COMMON NUVRO*NX*NY«NZ*NSEC* NSP1 

COMMON NC*NC1 *NC2*NC3 *NR*NR1* NR2*NR3*NR4 

COMMON NlttCyB 

COMMON ERR*TITLE*NFREQ 

GO T0(U*12*13«14),KD1 

11 RETURN 

12 KD2*5 
KD3-9 
BK*1. 

GO TO IS 

13 KD2*4 
KD3-9 
BK*1. 

GO TO 15 

14 KD2*5 
KD3*8 
BK*-1. 

15 N0UT*6 

GO TO (It 170*201) *NPIB 
1 DO 311 1*1,10 
DO 312 11*1*5 
BPL( I * 1 1 1*0* 

312 BND( I * II l*0« 

DO 313 11*1.7 

313 SAV( I * 1 1 )*0« 

311 CONTINUE 

BN0(2,1>*1. 

BND(KD2*2)*1. 

BND(6,3)*1. 

BN0(KD3.4)*l. 

BNDC 10*51*1* 

IF(NC-l) 103*103*100 

100 DO 101 I * 1 « NC 

101 SAV(I*I)*1. 

DO 102 I*2*NC 

102 SAV(I*1)*B(I-1> 

103 SAV(1*1)«1. 

C GO TO BRANCH IS NOT NEEOED 
105 00 112 M*NR1«NR2 
NN*0 

DO 110 N*NC1» NC2 
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NN*NN+1 

110 BPLfM*NCl ) *BNDf M* N) *SAV tNN* 1) +BPLfM*NCl) 

IFtNC-1) 120*120, 111 

111 NN*NC 1*1 

00 112 N*NN*NC2 

112 BPLfN,N)*BN0fM,N) 

120 NPIB*1 

GO TO 400 

170 IFfNC3.EQ.0IG0 TO 181 
00 180 N*1*NC3 
180 SAVtN+l*7)«BfN> 

181 X*i. 

GO TO 203 

201 X*BPLf 1*1) 

IF(NC3«EQ*0) GO TO 203 
00 202 N*1,NC3 

202 BfN)*SAVtN<>l*7) 

203 00 204 1*1, NR 

204 SAVf 1*11*0 
1*0 

210 IF(NRl-3) 211*212*213 

211 IFfNR2-6t 215*216*218 

212 IFfNR2-6> 220*220 ,222 

C EDGEWISE BENDING ONLY W*l. 

213 B t KD3 )*X/ f BNDf 7,4 l+BNDf 7,51 *Bfl)*BK) 

B( 10)*B( 1 )*Bf KD3) 

K04*K03— 1 
00 214 N*1,K04 

214 BtN)*0« 

IFIK04 •EQ*7) B( 9)*0 • 

NN*10 
GO TO 290 

C TORSION ONLY 

215 8t2)*X/BNDf 1*1) 

Bf 1>*0 

NN*2 

GO TO 290 

C TORSION FLAP 

216 VNORM*X/tBNDt 3,11+BNDf 3*2 )*B( l)+BN0f 3,3» *BI2 ) i 
SAV(2,1)*VN0RM 

SAVf K 02* 1 l*Bf 1 I ♦VNORM 
SAVf 6«l)*Bt2)*VN0RM 
00 217 N*1 *6 

217 Bf N)*SAVf N,1 ) 

NN*6 

GO TO 290 

218 VNORM*X/ ( BNDf 3 * 1 )+BND( 3* 2 ) *B ( 1 ) +BND ( 3 *3 ) *Bf 2 ) +BND( 3 *4 )*B f 3 1 ♦ 
1 BNO <3*51 *Bf4) ) 

SAVf 2 * 1 1 =VNORM 

SAVf KD2* 1 )*Bt 1 ) *VN0RM 
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SAV(6*1)*BI2)*VN0RM 
SAVIKD3*1)*BI3)4VN0RM 
SAVI 10*11 *BI 4)4 VNQRM 
00 219 N*1 • 10 

219 BIN)«SAVIN*1) 

NN*10 

GO TO 290 
C FLAP ONLY 

220 VNQRN*X/ IBNDI3*2)+BNDI3*3)*B(1)) 

BIKD2)*VN0RM 

BI6)*BI 11*VN0RM 
DO 221 N«l*4 

221 BIN)*0 
NN»6 

GO TO 290 

222 VN0RM*X/|BN0l3*2)*BN0t3f3)*Btl)+BNDI3*4)*BI2)+BND(3»5)*&t3) ) 
SAVI KD2 • 1 ) *VNORN 

SAVI6*1)*BI1)*VN0RM 
SAVIKD3* 1 )*VNORM 
SAVI 10* 1 )*BI 3 )*VNORM 

00 223 N*1 » 10 

223 BIN)*SAVIN*1) 

NN*10 

290 N*NN*'l 

1 FIN- 101 291*291 *293 

291 00 292 I*N*10 

292 BII)*0 

293 1*0 

WRITE I NQUT *920 ) I * I BIN) * N*l* 10) *XROOT 
BPLI 10*31*0 
BPLI 10*4) *XROOT 
920 FORMAT 1 13 *11 IE10 .3* IX) ) 

NPIB* I 

400 CONTINUE 
RETURN 
ENO 
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SUBROUTINE TEMAT II,T) 

TRANSFER MATRIX TRANSFERS ACROSS AN OFFSET 60+ ,E+5*LAS0+ 

A LUMPED MASSf AN ELASTIC FIELD* ANO A TWIST IN THE ELASTIC 
AXIS. SEE PAGE 13. 

SEE PAGE 15 FOR TERMS OF T MATRIX 

DIMENSION ELI 30 I * El XI 30)* El Y( 301 * El 2 C 30 1 *XINR( 30) .YINRI30) 
DIMENSION EMASt 30 )«PHII30)*EPSI30)*ELZI31)*CAPSII30) 

DIMENSION CCPSI30)*TI 10,101 *AVMASI31) 

DIMENSION BI10) 

DIMENSION TITLEI19)*ERR(10) 

COMMON CPSQ.OMSQ* XROOT* ZERMO 

COMMON EL*£IX«tt7*ElZ*XlNR«Yl NR»EMAS *PHI , EPS *ELZ,CAPSI «CCPS* EPSN 
COMMON CPOMG,OMEG,SIZER*THETC*ELC*AVMAS 
COMMON NGYRO*NX *NY*NZ*NSEC* NSP1 
COMMON NC «NC1 *NC2 *NC3 *NR*NR 1* NR2*NR3 *NR4 
COMMON NIBC*B 
COMMON ERR*TITLE* NFREQ 

IM1*I— l 
XI -XROOT 
HI»0. 

PTX»0. 

IFIIM1) 25*25*2 

2 DO 3 J-1*IM1 

HI IS DISTANCE IN Z DIRECTION BETWEEN PITCH AXIS AND PROJECTION OF 
MASS ON ROTOR PLANE. SEE PAGE 6B. 

Hl-Hl *ELZ I J^l l*CCPS I J ) 

3 X I-XI ♦ELI J) 

X-XI 
H-HI 

IF I I-NSPi ) 5,4,4 

4 SUMAS«.5*EMASINSEC) 

E PST* EPS I NSEC) 

EPS I NSEC I -EPSN 
GO TO 30 

C PTZ IS A FORCE IN THE Z DIRECTION. SEE PAGE 6B. 

5 PTZ*.5*lEMASIIMl)'tEMASII>)*IH+EPStIMl) *CCPSIIM1)) 

DO 10 J-I «NSEC 

X-X+ELIJ) 

IF(J-NSEC) 8,7*7 

7 SUMAS-.5*EMASI J) 

ELZIJ*1)«0. 

GO TO 9 

8 SUMAS«.5*IEMASI J)*EMASt J+l) ) 

H-H+ELZI J+l )*CCPSI J ) 

9 PTX»PTX*SUMAS*X 

PTZ-PTZ>SUMAS*I H+ I EPS I J >+ELZ( J*1 ) )*CCPS( J ) > 

10 CONTINUE 

PTX-CPSQ*PTX 

PTZ=CPSQ*PTZ 
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X-XI 

H»HI 

20 SUMAS*.5*IEMASI IMlI+EMASd)) 

30 AVMASI I )*SUMAS 

SNSi s SINI CAPS I ( INI) ) 

CSSI*CCPS(IM1) 

SCSI»SNSI*CSSl*CPSQ 
SNSQ*SNSI *SNSI*CPSQ 
CSSQ*CSSI*CSSI*CPSQ 
EPSEM-EPS4 IM1I*SUMAS 
H*H*EPS4IM1)*CSSI 
F6*SUMA$*CPSG*H 
F10»-F6*CSSI 
F6*F6*SNSI 

F2»YINRI IHil*SCSI>EPS( IM1 )*F6 

A5l—EPSEM*X*CPSQ 

A58»-XINR4 I Ml 1*SCSI“EPS4 I Ml )*F6 

A63»SUMAS *4 OMSQ+SNSQ) 

A67»“SUMAS *SCSI 
A6i*-EPS4 IH1I*A63 

A98*-XINRI IM1 )*SNSQ“Y INRI IM1) *OHSQ“EPSEM* 

1 EPS! IM1I*I0MSQ+CPSQI“EPSI IMil*F10 
A107* SUHAS *40MSQ+CSSQ) 

970 FORMAT I //I 5E20. 12 ) ) 

C MASS TERMS 
Td v l)«l. 

TI2*2)«1. 

T (2*3 )*-A61 

Tl 2*1 1»-0MSQ*XINR lIMll + YINRIIMll* ( C SSQ-SNSQI -EPS ( I Ml ) * 
1 € T € 2*3)+F10) 

T42*7)«EPSUM1)*A67 
IFII.GE.NSP1) GO TO 35 
TI3,3)*COS4PHId) ) 

T 4 3* 7 I »—S INI PHI 111) 

GO TO 36 

35 T43*3)«1.0 
T43*7)»0.0 

36 CONTINUE 

TI6«3)sTI 3*3I*A63-*T43*7 )*A67 
T46«l I*- EPS 4 IMil*TI6f31 
T 1 6*6 l»T I 3,3) 

T 16,7 )«TI 3,3I*A67+T 43*7 )*A107 

T46*10)*TI3*7I 

T4 7*3 I»-T 4 3*7) 

T 1 7*7 ) »T 4 3 *31 

TI10*3I«T 13* 3I*A67-T 1 3*7) *A63 
Tll0*ll»-EPS4IM1I*T410*3) 

T 1 10,6) *-T 4 3,71 

TI10,7)*T4 3*31 *A 107- T43*7I *A67 
T 4 10*10)*T (3*3) 
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49 IFII-NSP1I 50*450 *450 

25 HI-ELZtl) 

H*HI 

XI-XROOT 

X*XI 

C PTZ IS THE FORCE IN THE Z DIRECTION 
PTZ«.5*EMASI1)*H 
DO 29 J*1 *NSEC 
X*X*EL(J) 

I FI J-NSEC) 27*26*26 

26 SUMAS*.5*EMASIJ) 

GO TO 28 

27 SUMAS*.5*IEMASt J)+EMASIJ*1) I 
H*H«-ELZI J+1)*CCPSI J) 

28 PTX*PTX^ SUHAS*X 

P T Z*P T L* SUHAS* I H* IEPSI J l + ELZI J+l ))*CCPSIJ) I 

29 CONTINUE 
PTX«PTX*CPSQ 
PTZ*PTZ*CPSQ 
X*XI 

H-HI 

C 

31 SUMAS*.5*EMAStl> 

AVMASI1)*SUMAS 

F10»-SUHAS4CPSO*H 

A63*SUMAS*QMSQ 

A107*SUNAS4|0MSQ*CPSQI 

A51»0. 

A 58*0* 

A61*0* 

A67»0. 

A 98*0 • 

C MASS TERMS 
Tl 1*1 1*1 • 

T (2*21*1. 

T(3*3)*C0SIPH1C II) 

T(3*7)*SINIPHII I) ) 

T(6*3)«T<3*3>*A63 
T (6*6 )*T 13*3) 

TI6*7I*T 13*71 *A 107 
T I6«10)*T 13*7) 

TI7*3 l*-T I 3*7 1 
Tl 7*7 )*T (3*3) 

T I 10*3)*-T (3*7) *A63 
TI10,6I*-TI3,7) 

TI10*7)«TI3*3)*A107 
T I 10* 10)*T 13*31 
GO TO 49 

50 I FI NX) 100*100,150 
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100 T(l t 2l*EL< I>/Eixm 
TIl#3l*TIl#2l*T(2,3) 
T(1,1I*H1,1)+T(1»2)*T(2»1) 
Tll,7»*T(l#2>*Tf2,7» 


V AXIS TERMS 
150 IFINYI 160# 160* 180 
160 L*1 
165 EI-EIY(I) 

C SEE PAGE 6B FOR DEFINITION OF GAMMA 

1 TO GAMMA*EL I i l*SQRT t PTX/EI ) 
G2»GAMMA*6AMMA 
IF (G2~«00001 1172# 172# 171 

171 SOC-ITANHI GAMMA >1**2 
S«SQRT<SOC/<l.-SOCn /GAMMA 
C»SQRT(l./<l.-SOC)> 

GO TO 173 

172 S«l.+G2/6. 

C«i.*G2*.5 

173 C1«S*T(10#6) 

C2*S*T(10#10) 

C3«C*T<10,6| 

C4*C*T<10#10) 

IFIG2-. 000011 174# 174# 168 

174 S-I1.+G2*. 051/6. 

C«.5*G2/24. 

GO TO 169 

168 S»(S-1.I/G2 
C*(C-1.)/G2 

169 C5»S*TI 10,6) 

c6«s*mo,io» 

C7*C*T ( 10,6) 

C8«C*T< 10,101 

SL»ELU)/EI 

SL2-EL( I)*SL 

SL3»EUI»*SL2 

GO TO (175,178# 175#178)#L 

175 CL»SL*Cl 
CL2*SL2*C7 
CL3»SL3*C5 

GO TO ( 176 #230# 332,360) ,L 
C Y AXIS - SINE TERMS 

176 TC3,ll»*CL3*T(2,7) 

T ( 3 ,3 1*-CL3*A67+T 13,3) 

T(3,7I»-CL3*A107*T(3,71 

TC5,10J*-EL(I)*C1 

T (3,8)«~CL2*A98+T(5,10I 

T(3,9I»-CL2 

TI3,10»=-CL3 
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T(4il )*+CL2*T(2«7 ) 

T(4,3)»-CL24A67 

TU#7I»-CL2*A107 

7(4,8 »«-C3-CL*A98 

T(4»9)*-CL 

T(4,10)*-CL2 

T( 5*1 1*-T (5*10) *T (2,7) 

T(5,3)*T( 5,10)4 A6 7 

T(5,7)*T(5,10)*A107 

T ( 5,8 )*PTX*T (5,10 )—C3*A98 

T ( 5»9)*-C3 

178 CL»SL*C2 
CL2«SL2*C8 
CL3*SL3*C6 

GO TO (179, 210, 33 5, 357), L 
C Y AXIS - COSINE TERMS 

179 T ( 7»1 l*-CL3*T(2,7 ) 

T ( 7,3 l»CL3*A67*T( 7, 3) 

T(7,7)»T(7,7I«-CL34A107 

T(9,10)«EL(I)*C2 

T ( 7«8)*CL2*A98+T (9,10) 

T(7,9)»CL2 

T(7,10)»CL3 

T(8,3)*CL2*A67 

T ( 8,7)»CL2*A107 

T(8,8)«C4*CL*A98 

T(8,9)«CL 

T(8«10)*CL2 

T(9,3 )*T ( 9,10)*A67 

T(9,7)»T(9,10)*A107 

T(9,8)*C4*A98+PTX*T(9,10) 

T(9,9)*C4 

IF(IMl) 200,200,177 

177 T(8,1)*-T(8,3)*EPS(IM1) 
T(9,l)*-T(9»3)*EPS(IM1) 

C Y AXIS INFINITELY STIFF TERMS 

180 T(3,8)*-El(I)*T(10,6) 

T (4,8)»— T (10,6) 

T ( 5,1 )«-T (3,8) *T (2,7) 

T(5,3)*T( 3 ,81 *A67 
T(5*7)»T(3,8)*A107 
T(5,8)«PTX*T(3,8)-T(10,6)*A98 
T(5,9)*-T(10,6) 

T ( 5, 10)*T (3,8) 

T(7,8)*EL(I)*T(10,10) 

T(8,8)*T(3,3) 

T(9,3)»T(7,8)*A67 
T(9,7)»T(7,8)*A107 
T(9,8)*T( 10,10) *A98+PTX*T ( 7,8 ) 
T(9,9)=T< 10,10) 
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TI9*101«TI7,8I 
I Ft INI 1 200*200*181 
181 Tt9,ll*-EPStIMll*A67*Tt7,81 

Z AXIS TERMS 

200 IF4NZ1 205*205*250 
205 EI«EIZ4II 
L*2 

60 TO 170 
C Z AXIS - COSINE TERMS 

210 1(3*1 l*Ct2*A51-£L3*T 12*31+743*1) 
T43*31*T4 3*31 +CL3+A63 
T43*41«EL4 II*C2 
T 1 3*5 1*CL2 
T 4 3*6 l*CL3 

T 4 3*7 )»CL3*A67+Tt 3*7» 

T 4 3*81 a CL2*A58+T (3*8) 

T44,ll*CL2*A61+CL*A5i+TI4,l I 

TI4,3I»CL2*A63+T44,31 

T44,41»C4 

T(4,5I»CL 

T44,6)«CL2 

T 4 4*7 1»CL2*A67+T 14*71 
T4 4*8 1*CL*A58+T4 4*8) 

T4 5* 1 l«C4*A51-T 45*61*T(2*31 +TI5.1 

T 4 5*6 1»T 4 3*4) 
T45*31*A63*Tt5*6l+T45,3) 

T 4 5*41*PTX*T4 5*61 
T 4 5*5 )»C4 

T45*7I»T 4 5*61*A67+T 15*71 
T 1 5*8 1*C4*A58+T 15*8) 

60 TO 175 

C Z AXIS - SINE TERMS 

230 T 1 7, i l*CL2*A5 1+CL3*A6 1+Tt 7, 11 
T 4 7*31*CL3*A63+T 4 7*31 
T47*4I»EL4 I1*C1 
T I 7.5 1*CL2 
T47,61»CL3 

T47,7I*CL3*A67+T4 7*71 

T I 7*8 1 »CL2*A58+T 47*81 

T I 8*1 1*CL*A51+CL2*A61+T 4 8*11 

T 4 8,3 l=CL2*A63+T { 8*31 

T 4 8»4)*C3 

T 1 8*5 1*CL 

TI8*61«CL2 

T 4 8,7 1«CL2*A67+T I 8,71 

T 4 8*8 1=CL*A58+T (8*81 

T4 9*l 1*C3*A51+T (7,41*A61+T49» 11 

T49,5)=C3 
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T(9*6)»T4 7*4) 

T(9*3)«T(9*6)*A63+T(9,3) 

T (9*4)»PTX6T (9*6) 

T(9,7)«T(9,6)*A67+T(9,7) 

T (9*8 )»C3*A58*T ( 9*8 ) 

GO TO 280 

C l AXIS INFINITELY STIFF TERMS 
250 T(3,4)»EL(I)*T(10,10) 

T ( 4*4) *T( 10*10) 

T45,1)*T( 10*10)*( A51*EL(I )*A61)+T(5*1) 

T 4 5 ,3 )*T 4 3*4) ♦A63+T (5*3) 

T(5*4)*PTX*T(3«4) 

T(5,5)*T(10,10) 

T( 5*6 )*T (3*4) 

T45*7)*T4 3*4)* , A67aT (5*7) 
T(5*8)«T(10*10)«A58+T(5*8) 

T(7*4)«EL( l)*T(10,6) 

T ( 8*4)*T ( 10*6) 

T(9,1)»T(10,6)*< A51+EL(I)*A61)+T(9,1 ) 
T(9*3)«T(7*4)*A63+T(9*3) 

T 49*4)»PTX*T 17*4) 

T(9,5)-T( 10,6) 

T 4 9*6 )*T 17*4) 

T 19*7 )»T I 7*4)*A67+T (9,7) 

T49*8)»T4 10*6)*A5 8«-T 19*8) 

280 IF(NGYRO) 400*400,300 

GYROSCOPIC TERMS 
300 OMEGA*SQRT (OMSQ) 

A68*-2* *EMAS4 IM1 ) *CP0MG*EPS4 1 M1)*QMEGA 
A97*— A68*CSS I 
A68«A68*SNSI 

A28*-(XINR4 IM1 )*Y INRI IM1) )*CP0MG*SNSH-A68 
A91—A28 

T42*8)*A28 

T(6*8)»T( 10,10)*A68-T(10,6) *A97 
T410*8)«T410»6) *A68+T (10, 10 )*A97 
IF ( I -NSEC) 305*500*500 
305 IF(NX) 310*310*320 
310 T41,8)*4EL4 D/EIXf I))*A28 
320 IF (NY) 330*330*340 
330 L*3 

GO TO 165 

332 T43,i)*T43,l)-CL2*A91 
T43*3)»T( 3*3)-CL2*A68 
T(3,7)*T( 3,7)-CL2*A97 
T(3*8)*T( 3,8)-CL3*A97 
T(4,l)«T(4,l)-CL*A9i 
TC4,3)=T(4,3)-CL*A68 
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T 44*7 1*T 4 4*8l-CL4A97 

T 44*8 1*T4 4*8 l-CL2*A97 

T45*il*-C3*A91 

T(5*3)«- C3*A68 

T45*7)»-C3*A97 

T45«8I*-EL4I1*C1*A97+T45*8I 

GO TO 178 

335 T47*1I*T47«11*CL2*A91 
T47*31*T4 7»31+CL2*A68 
T47*71*T47,71*CL2*A97 
T4 7*8)«T4 7*8I«>CL3*A97 
T(8«l)«Ti8'll*CL*A9l 
T4 8*3 J*T 4 8«31+CL*A68 
T48.71«T<8.7M-CL*A97 
T48*81«T48*81+CL2*A97 
T49*1)*T49,11+C4*A91 
T49.3I-T49.3I*C4*A68 
T49*71«T49*7I+C4*A97 
T49*81*T49 V 81+EUI1*C2*A97 
GO TO 350 

340 T45*1I«T45*1I~T410*61*A91 
T45*31«T45*31-T4 10*61 *A68 
T45*7I»T45.71-T410*61-A97 
145*81 *T 4 5*81-EL4 I l*T4 10*6) *A97 
T49*11*T410*101*A91*T49,1) 

T(9*31»T410*10 1*A68+T 4 9* 31 
T<9,71-T49.7)*T410,10)*A97 
T (9,8)»T( 9*81 *Et.4 1»*T4 10*10 l*A97 

350 IF4NZ 1 355*355*380 

355 E I»EIZ4 1 1 
L*4 

GO TO 178 

357 T410*8)*T4 10* 81 ♦CL3*A68 
T44*81*T44»81*CL2*A68 
T45«81*T45«81*EL4 II*C2*A68 
GO TO 175 

360 T47*8l3T4 7*81 ♦CL3*A68 
T ( 8*8 )*CL2*A68+T 48*8) 

T49*8)«EL4II*C1*A68 
GO TO 400 

380 T45*8)«T45*81+EL4 I l*T 4 10. 10 1*A68 
T49,81=T49.8»+EL( I l*T4 10* 6) *A68 
ENO GYROSCOPIC TERMS 

400 IF4ELZ4 I I.EQ.01 GO TO 500 

MODIFY T MATRIX FOR RIGID OFFSET BY POST MULTIPLYING BY OFFSET MAT 
SEE PAGE 10 

410 C*PT Z*ELZ 4 I 1 
DO 420 J»l*lO 
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Tl J»1I*TI4*1I*C*T ( J*2 l-ELZI I l*T(J»3) 
TI4*6)*T( J*6)-ELZ< I )*T(J*2) 

420 T< J,8I«T( J,8>*C*T( J,9) 

GO TO 500 

REMAINOER OF NASS MATRIX AT LAST SECTION 
450 00 457 J*l*10 
457 Tt4*4)*l. 

EPSN*EPS4 NSEC ) 

EPS(NSEC) *EPST 
T(5*l )*A51 
T 1 5*8 I*A50 
T(9*8)*A98 
500 RETURN 
ENO 
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SUBROUTINE MATINV ! A, N1 ,8, Ml, DETERM, ID) 

DIMENSION A! 10. 10 1 • B( 10. 11 • INDEX! 10, 31 

EQUIVALENCE < IROW.JROW) .( ICOLUM, JCOLUM) , ( AMAX, T, SWAP ) 

M*1 

N»N1 

10 DETERM*1 • 

15 00 20 J*1 ,N 
20 INDEX! Jf 31*0* 

30 00 550 I»1 ,N 
40 AMAX*0. 

45 00 105 J»1,N 

IF! INOEX! J.3)-l)60, 105*60 
60 00 100 K*1.N 

IF ! INDEX!K,3)-1 ) 80.100.715 
80 IF ! AMAX-ABS! A! J.K) ) ) 85.100.100 
85 IROU-J 
90 IC0LUM«K 

AMAX-ABS! A! J.K) ) 

100 CONTINUE 
105 CONTINUE 

INDEX! ICOLUM, 3)* I NOEX ! ICOLUMt 3) +1 
260 I NOEX! I. 1 )*IROM 
270 1 NOEX ! I . 2 ) "ICOLUM 
130 IF ! I ROW* ICOLUM ) 140.310.140 
140 DETERM*~DETERM 
150 00 200 L*1,N 
160 SMAP*A! IRON. LI 
170 A! IR0N«L)*A! ICOLUM. L) 

200 A! ICOLUM. L)«SWAP 
IF CM) 310.310.210 
210 00 250 L»1,M 
220 SHAPES! IROM.L) 

230 8! I ROW. L 1*8! ICOLUM. L) 

250 8! ICOLUM, L)*SWAP 
310 PIVOT*A! ICOLUM, ICOLUM) 

OETERM*OETERM*PIVOT 
330 A! ICOLUM. ICOLUM)* 1* 

340 00 350 L*1,N 

350 A ! ICOLUM, L )* A! ICOLUM. L) /PI VOT 
355 IF CM) 380.380,360 
360 00 370 L«1,M 

370 B! ICOLUM, L)*BC ICOLUM, D/PIVOT 
380 DO 550 Ll»l,N 
390 IF ILl-ICOLUM) 40C, 550,400 
400 T*A!L1, ICOLUM) 

420 AILl. IC0LUM)=0. 

430 00 450 L»l,N 

450 A!L1,L)*A! LI , L) -A (I COLUM, L ) *T 
455 IF (M) 550,550,460 
460 00 500 L=l,M 
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500 BUltU*B<Ll«U-B(ICOLUM,U*T 
550 CONTINUE 
715 10*2 
740 RETURN 
ENO 


SUBROUTINE FINOH4 0MLST*0MSQ*DET*DT0LD*RNEW*RQLD*NEXP0*NC4* NFREQ*OM 
1INT*NTIME*NFIN0*NTRL«NFU,,NPDF) 

c 

DIMENSION RNEWI4) «R0LD(4) *OMINT( 15) *0M(4) ,0MRD(4) 

C 

900 FORMAT (8(5X*I5)) 

910 FORMATI 4X * 13 *3X* E20. 12*E15.5*15X*4E15.5) 

920 FORMATI 4X* I3*3X* E20. 12*6E15. 5) 

C 

N0UT»6 

IFINF IN0.GT.100) GO TO 130 
1F(NFLG.GT .01 GO TO 1 
OMLSS-1. 

NFLG*10 

DET1-0ET 

IF(NPOF.LE.O) GO TO 3 
00 2 I-l.NPOF 

2 OMINT 1 1 )*OMINT( I ) *OMINT II) 

3 NFREQ»NPDF*1 

1 IFINTIME.NE.O) GO TO 10 
DET*0ET1 
OMLST-1. 

OMSQ-l.OOl 
OMSQ*OMSQ* 1.001 
GO TO 6 

5 OMEGA » SORT I ABS( OMLSS) )*OMLSS/ABS( OMLSS ) 

IF INCA. NE.O. ANO.NFIND.LT. 1001 

1WRITE 16*920) NTIME. OMEGA. DET *YB* ( RNEW ( J ) .J-1.NC4) 

IFINC4.EQ.0. ANO.NFIND.LT. 1001 
1MRITE (6*920) NT) ME* OMEGA « DET *YB 
IFINC4.EQ.0. AND.NFIND .GE. 100) 

1WRITE (6*920) NTIME*OMEGA*DET 
IFINC4. NE.O. AND. NF1ND.GE. 100) 

1WRITE(6*910) NT I ME* OMEGA* DET* (RNEM( J ) , J-l ,NC4 ) 

6 NTRL*2 
NTIME«NTIME+1 
RETURN 

10 OMD»i. 

OMOLD-1. 

NF»NFREQ-1 
IFINF.EQ.O) GO TO 50 
DO 30 I»1*NF 

QMOLD*OMQLD*(OML$T-OMINT( I ) ) 

30 OMD at OMD*(OMSO”OMINT ( I ) ) 

0 MP*0 MD*Q MOL. D 
IF(OMP.NE.O) GO TO 50 
OMLST*OMSQ 
0MSQ*0MSQ+1 • 

GO TO 5 
50 OMO*i «/OMD 
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OMOLD-l./OMOLD 
IFINFINO.GT.lOO) GO TO 130 
60 YA*DET*OMO 

YB*OTOLO*OMOLO 

OMORsOMSO 

OMLSS-OMLST 

CALL INTERPIOMLSS tOMOR* YA*YB) 

IF (ABS(OMOR-OMLSS)/OMLSS.LT.*316228**NEXPO) GO TO 70 

OMSQ»QMQR 

OHLST-OMLSS 

GO TO 5 

TO NFlNO»NF INO^lOO 
0TSAV*0£T 
OMSAV»OMSQ 
OMNXT-OMOR 

IF (NC4.EQ.0) GO TO 130 
00 90 I»1*NC4 
VA*RNEWI I ) 

YB*R0LD4 I ft 

OMOR»OMSQ 

OMLSS-OMLST 

CALL INTERP(OMLSS*OMOR» YA»YB) 

OMIII-OMOR 

90 QMRDI I l-ABSIOMNXT-OMOR) 

NOOR-1 

0MQ-0MR0I1I 

00 110 I-1.NC4 

IF (OHQ.LT.OMRO(II) GO TO 110 
NOOR-I 
OMQ-OMROIII 
110 CONTINUE 
OMLST-OMSQ 
OMSQ-OMCNOORI 
WRITE 16*9301 NOOR 

930 FORMAT! /* 5X.43HSUBSEQUENT ITERATIONS USE REMAINDER NUMBER *13) 
GO TO 5 
130 OMOR-OMSQ 

IF INTIME.LT.2.AND.NC4.NE.0I GO TO 60 
OMLSS-OMLST 

IF INC4.NE.0I GO TO 135 
YA-OET 
YB-DTOLO 
GO TO 140 
135 YA-RNEWINOORI 
YB-ROLDINQOR) 

140 CALL INTERPIOMLSS*OMOR*YA*YBI 

1 FI (QMOR-QMLSS) /OMLSS »LE . • 1**NEXPQ) GO TO 150 
OMLST =OMSQ 

OMSQ*OMOR 
GO TO 5 
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150 NF I ND»NF I NO— 100 
ONINT(NFREQ)*OMOR 
NTRL*i 
NTIME*0 
RETURN 
END 



SUBROUTINE INTERP < OMLST , OMSQ,YA,YB> 
N0UT*6 

0M*0MSQ-YA*40MSQ-CMLST)/4 YA-YB) 

IF(ON) 1,2,2 

1 WRITE (NOUT,900) 
omlst»omsq 
OMSQ»QMSO*1.1 

900 FORMAT 4 30H ITERATION GIVES OM NEGATIVE 
GO TO 3 

2 OMLST-OMSO 
OMSQ*OM 

3 CONTINUE 
RETURN 
ENO 
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SUBROUTINE SHAPE < PRO* NT APE * NPCH, I BC • TMSL } 

DIMENSION EL(30)»E1X(30),EIY(30).EIZ(30) ,XINR(30) . YINRC30 ) 
DIMENSION EMAS(30)«PH1(30) «EPS(30),ELZ(31) .CAPS I (30) 

DIMENSION CCPS(30),AVMAS<31 ) 

DIMENSION PRDI10,5),T(10,10).Z( 10),ZSV< 10) ,S( 10) 

DIMENSION A(10«10),C( 10) 

DIMENSION ERR (10) , AS AVI 10* 5 )t TITLE! 19) 

COMMON CPSQ,OMSQ« XROOT.ZERMO 

COMMON EL.EIX.EIY.EIZ ,XINR, YINR.EMAS ,PH1 *EPS , ELZ»CAPSI ,CCPS , EPSN 

COMMON CP0MG,0MEG,S1ZER,THETC«ELC«AVMAS 

CCMMON NGYR0,NX,NY«NZ«NSEC,NSP1 

COMMON NO .NCI ,NC2 .NC3 .NR, NR 1, NR 2 , NR 3* NR A 

COMMON NIBC.B 

COMMON ERR«TITLE,NFREQ 

C SHAPE SUBROUTINE FOR PIN AT 0 AND LEAD-LAG HINGE AT POINT 1 
NQUT>6 
TMR-0. 

WRITE (NOUT .910) 

C INITIALIZE STATE VECTOR 

DO 10 11*1,10 
ZSVIII)*0, 

10 Z(II)*0. 

N*2 

IF(IBC.EQ.l) GO TO 401 
CALL BNDRVXIPRD,A,N, IBC ) 

GO TO 405 

401 CALL BNDRYKPRD, A,N) 

405 K*1 

OMEGA* SORT(OMSQ) 

16 CONTINUE 
GNMS*A( 10,5) 

X*A( 10,4) 

I*N 

DO 17 J*NR1,NR2 
Z( J)«B( J) 

17 ZSVU)*0. 

IF (K.EQ.l) GO TO 52 
IFITMR.GT.TMSL) GO TO 50 

IF(NPCH.NE.O) WRITE 14.916) OMEGA, ( T ITLE4 IKJ) . IKJ*1 , 16 ) .NFREQ 
GO TO 50 

COMPUTE AN OUTPUT STATE VECTOR AT EACH SECTION, AND COMPUTE 
GENERALIZED MASS, AND WRITE ON TAPE FOR EACH SECTION ASSOCIATED 
WITH A CONTRIBUTION TO THE GENERALIZED MASS 
35 DO 40 J»NRi,NR2 
Z(J)*ZSV( J) 

40 ZSV(J)*0. 

42 C(1)*AVMAS(I«-1)*EPS(I) 

C(2)*C( 1)*EPS< I ) 

GNMS*GNMS+(XINR(I )+C< 2) l*< Z <1 )*Z< 1)+Z(8)*Z(0) )-2.*C( 1)*Z(3)*Z( 1) 
1+ AVMAS ( lFl)*(Z(3)*Z(3)*Z(7)*ZC7) ) 
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IF4Z43I.EQ.O) GO TO 49 
IF 4K.EQ.2.0R.I.LT.NSEC) GO TO 49 
TMR*ABS4Z41)/Z43) I 
49 CONTINUE 
X»X*EL4I) 

WRITE 4NOUT*920)I*4Z4J)*J-1*10)*X 
50 IF4K.EQ. 1 .OR.NPCH.EQ.O.QR. I .EQ.O) GO TO 52 
51 IF 4TMR.LE.TMSL) WRITE44*930) 

1 Z43)*Z(7)*Z41)*Z48)*Z(4),Z42)*Z45)*Z46) *Z49) *Z410) 

52 I«I*l 

IF 4 1 -NSEC 1 150* 15Cf 500 

150 00 151 11*1*10 
00 151 IX«1*10 

151 T4 1 1* IX)*0. 

CALL TEMAT4 I *TI 

IF 4 I -NSEC) 170* 16C«500 
C ACCOUNT FOR THIS MASS AT OUTBOARD END 

160 00 161 M*NR1,NR2 
00 161 J*NR1,NR2 

161 A4M*J)*T4M,J) 

00 1611 11 * 1,10 

00 1611 IX«1*10 

1611 T4 1 1* IX)*0. 

I*NSP1 

CALL TEMAT 4l*T) 

I»NSEC 

00 163 M*NR1«NR2 
00 162 J*NR1«NR2 
C4 J)*0. 

00 162 N*NR1»NR2 

162 C4J)*C4 J)+T4M*NI*A4N*J) 

00 163 J*NR1*NR2 

163 T4M*J)*C4J) 

C MULTIPLY STATE VECTOR 41) BY T TO OBTAIN STATE VECTOR 41+1). 

170 00 180 J*NR1*NR2 
00 180 N*NR1*NR2 
180 ZSV4J)*ZSV4J)*T4J*N)*Z4N) 

GO TO 35 

C WRITE NATURAL FREQUENCY IN RAO/SEC ANO X CAP-OMEGA 

500 X*OMEGA/CPOMG 
GNMS* SORT ( GNMS ) 

WRITE 4N0UT *915) OMEGA* X*GN MS 
GNMS* 1. /GNMS 

IF4TMR.GT.TMSL.AND.K.EQ.2.AND.NPCH.EQ.1) WRITE4NOUT*917) 

C COMPUTE MODE SHAPE FOR GENERALIZED MASS IF NEEDED 

501 GO TO 4 502*510) *K 

502 K*2 
X«XROOT 

506 A! 1 « 1 ) =GNMS 

C COMPUTE GENERALIZED MODE SHAPE 


C COMPUTE GENERALIZED MASS MODE SHAPES THROUGH NLBC IN SND8Y IF 

C NIBC GT. 2 

N=3 

IP (IBC. EQ. 1) GO TO 411 
CALL BFDRYX (PRD, A, N, IBC) 

G O TO 1 6 

411 CALL BNDKY 1 (PRD, A,N) 

GO TO 1u 
510 CONTINUE 

910 FORMAT (1 (/) , 1 


1 2 1 H I 

ANGLE 

TORQUE 

DISPL 

A LZ 

MOM ENT 

SH 

2EaR 

DISPl 

ANGLE 

MOMENT 

3 HfiAEt 

RADIAL/ 

1 

3 2 i ri 

PHI 

T 

'v 

THETA 

M 2 

— 

4 V Y 

-W 

PS I 

H Y 

V Z 

COORD. ,/) 



920 FORMAT (£3, 1 1 (£10. 3, 12) ) 

930 FORMAT (5 ( 31 4. 7 , 1 X) ,/, 5£15. 7) 

915 FORMAT (/ 9 X , 1 9 H NATURAL FREQUENCY =F1 0. 4 , 2X , 1 4 HI AD IANS/SECON D, 

1 3H = , F I 0. 4, 2 OH X ROTATIONAL SPEED. 3 X ,20HGEN£RALIZED MASS 

2 = , F 1 0 . 4 //) 

916 FORMAT (E 12. 5, 16A4, 12) 

917 FORMAT (//3X,107H ** NOTE ** THE ABOVE MODE SHAPE IS NOT INCLUDED 
1 1N PUNCHED OUTPUT, PHI (TIP) .GT. (V (TIP) *TMSL, TORSION LIMIT) 

RETURN 

END 
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APPENDIX D 

TEST CASES AND ANSWER LISTINGS 


This appendix contains, in order: 

(1) A listing of computer program input which was used to 
generate various test cases at NASA Langley , 

(2) Answer listings, chosen from the printed output list- 
ing produced by the above input, and labeled by a 
Test Case number, and 

(3) Answer listings, including all the punched cards in 
the same order as produced by the above input (part 1) . 

The portion of the input which corresponds to the enclosed 
printed output may be identified by comparing blade model and pro- 
gram control parameters. In order to facilitate this comparison, 
the following observations and table should be useful. Each model 
has a title card and cards read by the NAMELIST option. There were 
19 models in the test case run. If the models are numbered sequen- 
tially from 1 to 19, the following table gives the model number 
which corresponds to the various test cases as chosen from the 
printed output. 


Identification numbers 


Test Case Model 

1. a 3 

l.b 4 

1. c 5 

2. a 7 

2. b 10 

3. a 15 

3 . b 15 

4. 16 


The first three modes were punched for model number 2 and the 
last five modes were punched for model number 14. 
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APPENDIX E 


PROGRAM CHANGES FOR OPERATION ON AN IBM 360/65 


The program as listed in APPENDIX C is not operational on a 
computer which uses less than 14 digit floating point numbers. 
Changes to the program listed in APPENDIX C were made, and resulted 
in satisfactory operation for two test cases. Numerical results 
did not agree exactly with CDC 6600 resulting and appeared to be 
more accurate, probably due to the use of 16 digit numbers rather 
than 14 digit numbers. 

The program was compiled on an IBM 360/65 by FORTRAN IV G 
LEVEL 20 during March and April, 1972, at the University of 
Rochester, in Rochester, New York. The changes necessary for 
proper operation with this compiler and computer are given as 
follows : 


1. "Comment" or remove the first two cards in the main 
program. 

2. "Comment" or remove the CALL DAYTIM (DATE) or replace it 
by an equivalent call to obtain the calendar date. If 
this call is replaced by an "equivalent" call, the 
format statement where DATE is written may need to be 
revised. This is number 898 in the main program. 

3. Add the following cards to all subroutines: 

(a) before existing non-executable statements, 

IMPLICIT REAL* 8 (A-H , O-Z) 

(b) immediately following existing non-executable 
statements, ABS (X) = DABS (X) 

TANH(X) = DTANH(X) 

SIN (X) = DSIN(X) 

COS (X) = DCOS (X) 

SQRT(X) = DSQRT(X) 

4. (a) In subroutine MATINV, after the cards specified in 

3(b) add the statement: CALL ERRSET (207,256,10) 

(b) In subroutine MATINV, Just before the RETURN and 
END cards, add the statement: 

740 CALL ERRSET (207,1,1) 

5. In subroutine MATINV, before FORTRAN statement number 330, 

add the statement: IF (PIVOT . EO. 0 . ) GO TO 740 
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6. In subroutine FINDM, replace the FORTRAN statement 
immediately after that numbered 30, i.e., OMP=OMD*OMOLD , 
by: IF ( (OMD.EQ. 0. ) .OR. (OMOLD.EQ. 0. ) ) OMP=0. 

7. The data must be punched on a compatible keypunch. 

8 . The data in NAMELIST presently begins with $ and ends with 
$ symbols, but the IBM 360/65 requires that it begins with 
$ and ends with $END symbols. 

For operation with the University of Rochester's WATFOR, 

Version 1, Level 1, January 1970 compiler, the cards specified in 
4(a) and 4(b) should be replaced by the following: 

for 4(a): CALL TRAPS (0,0,1000) 

for 4(b): 740 CALL TRAPS (0,0,0) 

Items 4 and 5 are required to allow execution with underflow 
errors which occur when MATINV is called to evaluate a determinant. 
The underflows occurred when the subroutine was evaluating off- 
diagonal numbers for use in solving a set of simultaneous equations. 
These off-diagonal elements are not used in defining the deter- 
minant. Thus, use of ERRSET results in no numerical errors, but 
does allow program execution to continue when such problems occur. 
Item 6 is required to avoid an overflow in defining OMP. These 
changes were required on the IBM 360/65 because overflows and 
underflows result for numbers with smaller magnitude exponents than 
on the CDC 6600 or CDC 6400. The use of ERRSET or TRAPS is compiler- 
dependent. 
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